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Modeling  shapes,  be  it  in  2D  or  3D,  is  a  fundamental  constitutent  of  computer 
vision  and  computer  graphics.  Shape  modeling  techniques  have  been  widely  used  in 
shape  synthesis,  shape  reconstruction,  recognition,  motion  tracking,  and  so  on.  In  this 
dissertation,  we  introduce  a  compact  and  versatile  geometric  shape  modeling  scheme 
that  can  model  a  large  class  of  shapes  and  is  amenable  to  stable  and  efficient  numerical 
implementations.  Geometric  models  are  traditionally  well  suited  for  representing 
global  shapes  but  not  the  local  detail.  However,  in  this  thesis  we  propose  a  powerful 
geometric  shape  modeling  scheme  which  allows  for  the  representation  of  global  shapes 
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with  local  detail  and  permits  model  shaping  as  well  as  topological  changes  via  physics- 
based  control.  The  proposed  geometric  models  are  a  blend  of  geometric  and  physics- 
based  models  and  are  in  spirit  "similar"  to  the  now  popular  deformable  superquadric 
models  but  differ  from  them  considerably  in  the  expressiveness  and  numerical  stability 
leading  to  greater  applicability.  .    .  , 

The  proposed  modeling  scheme  consists  of  representing  shapes  by  pedal  curves 
and  surfaces — pedal  curves/surfaces  are  the  loci  of  the  foot  of  perpendiculars  to  the 
tangents  of  a  fixed  curve/surface  from  a  fixed  point  called  the  pedal  point.  By  varying 
the  location  of  the  pedal  point,  one  can  synthesize  a  large  class  of  shapes  which  exhibit 
both  local  and  global  deformations.  We  introduce  physics-based  control  for  shaping 
these  geometric  models  by  letting  the  control  point  vary  and  use  a  dynamic  spline, 
that  is,  a  snake,  to  represent  the  position  of  this  varying  pedal  point.  The  model 
dubbed  as  a  "snake  pedal"  allows  for  interactive  manipulation  via  forces  applied  to 
the  snake.  We  extend  the  model  to  automatically  cope  with  topological  changes  by 
introducing  the  concept  of  hybrid  geometric  active  contour /surface  models  wherein 
the  traditional  snake  in  the  snake  pedal  is  replaced  with  a  geometric  active  contour 
which  is  realized  via  a  level-set  implementation. 

Efficient  numerical  algorithms  for  shape  recovery  from  image  data  using  the  pro- 
posed snake  pedal  model  are  presented.  These  algorithms  involve  novel  mathematical 
derivations  that  lead  to  efficient  numerical  solutions  of  the  model  fitting  problem.  We 
demonstrate  the  applicability  of  this  modeling  scheme  via  shape  synthesis  examples 
and  shape  estimation  examples  from  image  data. 
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CHAPTER  1 
INTRODUCTION 

Shape  modeling  is  a  fundamental  constituent  of  computer  vision  and  computer 
graphics.  In  computer  graphics,  it  can  be  used  to  synthesize  shapes,  object  mo- 
tion and  object  interaction  for  the  purpose  of  computer  aided  design  and  computer 
animation.  In  computer  vision,  it  is  very  crucial  in  2D/3D  shape  representation, 
shape  reconstruction,  quantitative  model  extraction  from  image  data  for  analysis, 
visualization,  recognition  and  motion  tracking.  Recently,  shape  modeling  techniques 
have  been  widely  used  in  medical  image  analysis  applications.  Two-dimensional 
and  three-dimensional  shape  models  have  been  used  in  X-ray,  computer  tomography 
(CT),  magnetic  resonance  (MR),  and  ultrasound  images  to  segment,  visualize,  track, 
and  quantify  a  variety  of  anatomic  structures  including  brain,  heart,  cerebral,  lungs, 
liver  and  skull. 

Two  primary  tasks  of  importance  in  computer  vision  are  shape  recognition  and 
shape  recovery  from  image  data.  The  shape  recovery  task  imposes  the  requirement 
that  the  shape  model  be  capable  of  representing  fine  detail  and  hence  requires  that  the 
model  have  a  large  number  of  degrees  of  freedom.  Shape  recognition,  on  the  other 
hand,  requires  that  the  model  representation  be  compact  in  terms  of  the  number 
of  parameters  so  as  to  achieve  easy  storage  of  a  large  database  of  models  and  easy 
matching  for  retrieval.  These  two  requirements  are  potentially  conflicting.  Therefore, 
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designing  shape  models  that  can  possibly  "satisfy"  these  requirements  has  been  the 
focus  of  many  researchers  in  the  recent  years  in  the  computer  vision  community 
[3,  11,  14,  54,  78]. 

1.1    Problem  Statement 

There  are  numerous  shape  modeUng  techniques  in  the  computer  vision  and 
computer  graphics  literature.  Most  modeling  schemes  fall  into  two  major  categories, 
namely,  active  models  and  passive  models.  Active  models  or  deformable,  physics- 
based  models  are  distributed  parameter  models.  Typical  examples  are  generalized 
spline  models  such  as  snakes  [31].  Active  models  have  broad  geometric  coverage 
and  hence  are  good  for  shape  recovery  and  reconstruction.  On  the  other  hand,  they 
are  not  well  suited  for  recognition  since  they  require  a  large  number  of  parameters 
for  their  representation.  Passive  models  or  geometric  models  are  lumped  parameter 
models.  Typical  examples  are  spheres,  cylinders  and  prisms  [4,  67].  Passive  models 
can  represent  shapes  with  very  few  parameters  and  they  are  good  for  recognition 
but  not  for  representing  detailed  shapes.  In  order  to  fill  in  the  gap  between  the 
requirements  of  reconstruction  and  recognition,  a  class  of  hybrid  models  have  been 
proposed.  For  example,  Terzopoulos  and  Metaxas  [78]  proposed  a  hybrid  model 
which  combines  the  descriptive  power  of  geometric  and  physics-based  models  via 
superposition  of  a  spline  on  a  conventional  superellipsoid.  Hybrid  models  are  a  com- 
promise between  the  two  extremes  and  occupy  both  ends  of  the  spectrum  in  shape 
description.  One  end  of  this  spectrum  is  occupied  by  lumped  parameter  models  and 
the  other  by  distributed  parameter  models.  However,  most  hybrid  models  introduce 
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a  lot  of  nonlinearities  into  the  model  representation  to  describe  global  deformations 
such  as  bending,  tapering,  and  twisting.  These  nonlinearities  increase  the  complexity 
in  the  model  representation  and  affect  the  numerical  stability  of  the  algorithms  used 
in  fitting  the  model  to  image  data.  Moreover,  the  hybrid  models  can  not  handle 
shapes  of  unknown  topologies  without  introducing  additional  nonlinearities  into  the 
modeling  via  blending  schemes  [14]  or  other  special  handling  schemes. 

In  this  thesis,  we  develop  a  compact  and  versatile  shape  modeling  scheme  with 
physics-based  control  for  shape  design  and  shape  recovery  from  image  data.  This 
modeling  technique  allows  for  the  representation  of  global  shapes  with  local  detail 
and  allows  model  shaping  via  physics-based  control.  The  local/global  descriptive 
power  of  the  model  satisfies  the  conflicting  requirements  of  both  shape  reconstruction 
and  shape  recognition.  In  addition,  topological  changes  can  be  handled  naturally  and 
automatically  via  a  level-set  implementation  in  the  proposed  model. 

In  the  remainder  of  this  chapter,  we  will  first  review  related  shape  modeling 
schemes  in  literature,  with  an  emphasis  on  the  advantages  and  disadvantages  of  each 
scheme,  and  then  present  the  overview  of  our  model. 

1.2    Literature  Review 

There  are  numerous  shape  modeling  schemes  in  computer  vision  and  computer 
graphics  literature. 
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1.2.1    Geometric  fc  Deformable  Models  with  Fixed  Topology 

Motived  by  the  generalized  cylinder  idea,  Kass  et  al.  [31]  proposed  a  deformable 
cylinder  model  constructed  from  generalized  splines.  They  developed  force  field  tech- 
niques for  fitting  their  model  to  monocular,  binocular,  and  dynamic  image  data. 
The  distributive  nature  of  this  deformable  model  enhances  its  descriptive  power  and 
allows  the  representation  of  natural  objects  with  asymmetries  and  fine  detail.  How- 
ever, the  generalized  spline  representation  of  the  model  does  not  explicitly  provide  a 
description  of  object  shapes  in  terms  of  a  few  parameters. 

Terzopoulos  and  Metaxas  [78]  developed  a  hybrid  dynamic  model — a  deformable 
superquadric  model  which  combines  the  descriptive  power  of  geometric  and  physics- 
based  models  via  superposition  of  a  membrane  spline  on  a  core  superquadric.  Vemuri 
and  Radisavljecvic  [81]  developed  a  multiresolution,  hybrid  modeling  scheme  which 
represents  the  displacement  function  of  a  deformable  superquadric  model  in  an  or- 
thonormal  wavelet  basis.  This  multi-resolution  basis  provides  the  model  with  the 
ability  to  continuously  transform  from  local  to  global  shape  deformations  thereby 
allowing  a  continuum  of  shape  models  to  be  created  and  to  be  presented  by  relatively 
fewer  parameters. 

Pentland  and  Williams  [54]  introduced  a  novel  shape  modeling  technique  based 
on  the  modal  representation,  which  can  be  used  to  represent  global  shapes  with  local 
detail.  The  number  of  parameters  required  for  this  representation  depends  on  the 
level  of  detail  needed. 
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Cohen  and  Cohen  [11]  and  Han  et  al.  [26]  proposed  a  hybrid  hyperquadric 
model  that  adds  an  exponential  term  to  a  hyperquadric  primitive.  The  addition  of 
the  exponential  term  to  the  hyperquadric  equation  allows  for  local  control  of  the 
shape  generated.  An  arbitrary  number  of  concavities  can  be  generated  with  this 
modeling  scheme.  This  modeling  scheme  differs  from  most  of  the  schemes  discussed 
above  in  that  this  scheme  is  basically  a  geometric  model  but  it  allows  for  global  as 
well  as  local  deformations.  The  model  representation  is  quite  compact  in  terms  of 
the  number  of  parameters  needed  to  describe  the  shape.  However,  the  computation 
of  local  concavities  is  nontrivial  and  human  interaction  is  needed  in  the  model  fitting 
process. 

Another  modeling  scheme  described  by  several  authors  in  [3,  27,  30,  59]  involves 
superimposing  a  global  volumetric  deformation — the  free  form  deformation (FFD) — 
on  the  core  superquadric.  Both  global  and  local  deformations  can  be  represented  by 
this  modeling  scheme,  however,  it  is  difficult  to  describe  the  shape  of  an  object  with 
complex  deformations  and  fine  detail  with  this  modeling  scheme. 

Recently,  O'Donnell  et  al.  [50]  developed  a  novel  solid  shape  modeling  scheme. 
This  model  includes  built-in  offsets  from  a  core  base  model.  Local  deformation  over 
the  scaled-offset  model  is  allowed  via  displacement  vectors  attached  to  the  default 
model.  Model  fitting  to  data  is  achieved  in  a  dynamic  framework  as  in  Terzopoulos 
and  Metaxas  [78].  This  modeling  scheme  is  well  suited  for  tasks  such  as  cardiac 
motion  recovery  from  tagged  MR  images. 
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1.2.2    Topologically  Adaptable  Models 

Surfaces  of  arbitrary  topology  were  modeled  using  particle  systems  to  model 
surfaces  in  Szeliski  and  Tonnesen  [72].  Particles  can  be  added  and  deleted  dynamically 
to  enlarge  and  trim  the  surface. 

An  interesting  and  topologically  adaptable  deformable  model  was  developed  by 
DeCarlo  and  Metaxas  [14]  via  the  use  of  parameterized  blending  functions.  Their 
scheme  can  determine  the  number  and  locations  of  the  holes  in  an  object  with  un- 
known topology  and  accordingly  change  the  topology  by  use  of  appropriate  error  of 
fit  criteria. 

Numerous  techniques  based  on  the  concept  of  level-sets  have  been  reported  in 
the  literature  [9,  10,  46,  58,  82].  These  techniques  allow  a  deformable  contour  not  only 
to  represent  long  tube-like  shapes  or  shapes  with  bifurcations,  but  also  to  dynamically 
sense  and  change  its  topology.  Another  modeling  scheme  called  "topologically  adapt- 
able snakes"  was  introduced  by  Mclnerney  and  Terzopoulos  [48].  This  novel  version 
of  snakes  allows  for  automatic  change  in  topology  during  shape  recovery.  A  common 
problem  associated  with  these  models  is  that  they  do  not  possess  a  mechanism  to 
characterize  the  global/core  shape  of  an  object. 

1.3  Contributions 

In  this  thesis,  we  propose  a  powerful  geometric  shape  modeling  scheme  which 
permits  the  representation  of  global  shapes  with  local  detail  using  only  a  "semi- 
global"  characterization  and  relatively  simple  and  efficient  numerical  methods.  The 


proposed  modeling  scheme  consists  of  representing  shapes  by  pedal  curves  and  surfaces- 
pedal  curves/surfaces  are  the  loci  of  the  foot  of  perpendiculars  to  the  tangents  of  a 
fixed  curve/surface  from  a  fixed  point  called  the  pedal  point.  By  varying  the  location 
of  the  pedal  point,  one  can  synthesize  a  large  class  of  shapes  which  exhibit  both  local 
and  global  deformations.  We  introduce  physics-based  control  for  shaping  these  geo- 
metric models  by  letting  the  pedal  point  vary  and  using  a  dynamic  spline,  namely, 
the  snake,  to  represent  the  position  of  this  varying  pedal  point.  This  process  of  gen- 
erating models  is  referred  to  as  the  "pedaling  operation."  The  model,  dubbed  as  a 
"snake  pedal,"  allows  for  interactive  manipulation  via  forces  applied  to  the  snake.  In 
this  model,  automatic  changes  in  topology  can  be  realized  by  embedding  the  snake 
in  a  level-set  based  evolution.  •  •  \ 

The  key  contributions  of  this  thesis  are  summarized  below  in  the  form  of  salient 
features  of  the  developed  modeling  scheme,  the  snake  pedal,  and  the  associated  nu- 
merical algorithms:  ■     .         ;  ' 

1.  The  snake  pedal  model  is  inherently  geometric  but  can  exhibit  both  local  and 
global  deformations  as  would  a  physics-based  model. 

2.  The  model  can  exhibit  large  global  deformations  such  as  bending,  twisting,  and 
tapering  without  explicitly  incorporating  parameters  for  the  same. 

3.  The  modeling  scheme  allows  (for  the  first  time)  the  concept  of  a  global  shape 
to  be  incorporated  into  the  curve/surface  evolution  framework. 
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4.  The  model  permits  a  level-set  implementation  and  thereby  seamlessly  allows 
for  topological  changes  in  shape  during  evolution. 

5.  Fast  numerical  algorithms  for  estimating  the  model  fits  to  2D/3D  data  sets. 

6.  A  hierarchical  neural  network  based  learning  scheme  for  incorporating  prior 
knowledge  about  shapes  of  interest  from  the  application  domain. 

We  elaborate  these  contributions  in  the  following  sections. 
1.3.1    Compact  Representation  of  the  Model 

The  snake  pedal  model  is  inherently  geometric  but  can  exhibit  both  local  and 
global  deformations  as  would  a  physics-based  model.  Our  modeling  scheme  is  some- 
what similar  to  the  hybrid  modeling  schemes  mentioned  in  Section  1.2.1  in  that 
both  techniques  combine  the  descriptive  power  of  physics-based  models  and  geomet- 
ric models.  The  difference  lies  in  that  the  hybrid  models  discussed  earlier  have  to 
explicitly  include  the  global  parameters  to  incorporate  global  deformations  such  as 
bending,  tapering,  and  twisting  to  achieve  the  required  deformations,  which  causes 
additional  nonlinearities  in  their  representation,  while  our  model  does  not.  This  com- 
pact representation  of  global  deformations  by  our  model  can  avoid  the  introduction  of 
additional  nonlinearities  into  the  model.  Therefore,  numerical  efficiency  and  stability 
can  be  largely  improved  in  the  model  fitting  to  the  image  data. 

This  compact  representation  comes  from  the  mechanism  of  creating  shapes  in 
the  snake  pedal  model.  The  pedaling  operation  is  a  nonlinear  operation.  When 
applying  the  pedaling  operation  to  geometric  primitives,  one  can  produce  global 
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deformations  of  the  geometric  primitives  by  performing  the  nonlinear  transformation. 
This  nonlinear  operation  is  only  applied  in  the  model  building  phase,  and  it  does  not 
affect  the  model  fitting  phase,  as  we  will  discuss  in  Chapter  5.  Hence  no  additional 
nonlinearities  are  added  to  the  model  fitting  process. 

1.3.2  Automatic  Topological  Change  Handling  Scheme 

The  proposed  snake  pedal  model  allows  for  the  estimation  of  shapes  of  arbitrary 
topology  without  resorting  to  the  blending  or  other  special  operations  as  required  in 
Szeliski  and  Tonnesen  [72].  The  topology  changes  are  achieved  naturally  by  embed- 
ding the  modeling  scheme  into  a  level-set  framework  [46].  We  introduce  the  novel 
concept  of  a  global/core  model  into  the  now  popular  partial  differential  equation 
(PDE)  based  curve/surface  evolution  framework.  This  global/core  model  allows  us 
to  express  any  shape  as  a  combination  of  a  global  shape  such  as  ellipsoid  and  sphere, 
and  a  variable  offset  with  respect  to  this  global  shape  and  an  evolving  curve/surface. 
Because  of  the  presence  of  the  global  shape,  the  snake  pedal  model  can  be  regarded 
as  a  unique  hybrid  geometric  active  model. 

1.3.3  Fast  Numerical  Methods 

Novel  and  fast  numerical  methods  are  proposed  to  achieve  stable  local  optimal 
solutions  in  the  model  fitting  process.  These  numerical  methods  are  not  a  regurgi- 
tation of  techniques  that  are  described  in  a  standard  numerical  analysis  text  book, 
or  for  that  matter,  the  book  on  numerical  recipes  in  C  [84].  In  all  the  methods 
that  are  discussed  here,  novel  derivations  are  presented  which  lead  to  development 
of  software  modules  that  can  either  be  easily  embedded  into  or  used  with  standard 
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numerical  software  such  as  the  conjugate  gradient  (CG)  code  and  the  alternating 
direction  implicit  (ADI)  code. 

1.4    Thesis  Organization 

The  remainder  of  the  thesis  is  organized  as  follows: 

Chapter  2  contains  a  review  of  the  basic  mathematical  foundations  of  classical 
active  models,  especially  the  snake  formulations  in  order  to  illustrate  the  physics- 
based  control  that  is  also  pertinent  to  the  snake  pedal  model. 

In  Chapter  3,  we  present  the  definitions  of  pedal  curves  and  surfaces  and  show 
examples  of  the  same  in  2D  and  3D,  followed  by  the  development  of  the  proposed 
model,  namely,  the  "snake  pedal"  with  a  demonstration  of  the  expressive  power  of 
the  modeling  scheme  via  illustrative  examples  of  shape  synthesis.  In  addition,  we 
also  discuss  the  ways  of  automatically  handling  change  of  topology. 

In  Chapter  4,  we  cast  the  model  fitting  process  in  a  probabilistic  framework, 
and  present  the  advantages  of  the  probabilistic  framework  over  the  deterministic 
approach.  We  will  devise  a  novel  supervising  learning  scheme  to  incorporate  the 
prior  knowledge  of  object  shapes  into  the  shape  estimation  from  image  data. 

In  Chapter  5,  we  present  efficient  and  stable  numerical  algorithms  for  fitting  the 
snake  pedals  to  real  image  data.  The  model  fitting  process  can  be  posed  as  a  nonlin- 
ear optimization  problem  and  can  be  dissected  into  global  and  local  shape  parameter 
estimation.  We  present  efficient  numerical  solutions  for  both  tasks,  and  then  com- 
bine the  global  and  local  procedures  appropriately  to  solve  the  overall  optimization 
problem. 
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In  Chapter  6,  we  extend  the  snake  pedal  model  to  cope  with  topological  changes 
by  introducing  the  concept  of  hybrid  geometric  active  models.  We  first  review  the 
traditional  geometric  active  curve  evolution  theory  and  its  level-set  implementation. 
Following  this,  we  present  our  hybrid  geometric  active  model  evolution  and  its  level- 
set  implementation.  The  hybrid  geometric  active  contour  evolution  involves  the  evo- 
lution of  both  snake  and  snake  pedal  curves/surfaces.  However,  by  employing  the 
insightful  relationship  between  the  two  curve/surface  evolutions,  we  solve  these  evolu- 
tions simultaneously  by  only  solving  for  one  of  them  and  thereby  dramatically  reduce 
the  computational  and  storage  requirements  when  applied  to  shape  recovery  tasks. 

In  Chapter  7,  we  present  the  numerical  techniques  used  for  computing  a  stable 
solution  to  the  hybrid  geometric  active  model  evolution  equation.  The  equation 
of  motion  for  the  model  can  be  posed  as  a  special  type  of  Hamilton-Jacobi  partial 
differential  equation  (PDE).  Based  on  the  traditional  level-set  approach,  the  motion 
equation  of  the  our  model  can  be  efficiently  solved  by  employing  entropy-satisfying 
upwind  finite  difference  plus  minmod  finite  difference  schemes. 

In  Chapter  8,  we  present  the  experimental  results  and  the  applications  of  the 
model  in  shape  modeling  and  recovery  and  conclude  in  Chapter  9. 

Appendix  A  gives  the  detail  derivation  of  the  stiffness  matrix,  which  corresponds 
to  the  internal  smoothness  energy  imposed  on  the  snake  in  the  snake  pedal  model. 
Appendix  B  contains  the  derivation  of  the  relationship  between  the  curvature  of  the 
snake  pedal  and  that  of  the  snake  at  a  specific  point. 


CHAPTER  2 

MATHEMATICAL  FOUNDATIONS  OF  ACTIVE  MODELS 


In  the  proposed  snake  pedal  model,  we  use  an  active  curve/surface  to  control 
the  model  shapes  and  deformations,  this  active  curve/surface  is  an  integrated  and 
non-separable  part  of  the  modeling  scheme.  In  this  chapter,  we  will  review  the 
mathematical  foundations  of  active  models  in  order  to  obtain  a  better  understanding 
of  the  behavior  of  the  controlling  curves/surfaces  in  the  snake  pedal  model. 

Active  or  deformable  curve  and  surface  models  gained  popularity  after  they  were 
proposed  for  use  in  computer  vision  [79]  and  computer  graphics  [77]  in  the  mid  1980s. 
In  computer  vision,  deformable  curve  and  surface  models  were  proposed  as  solutions 
to  ill-posed  inverse  problems  such  as  edge  detection  and  surface  reconstruction  [55, 
75].  The  theory  of  continuous  deformable  models  was  introduced  by  Terzopoulos 
in  a  Lagrangian  dynamic  setting  [76],  based  on  deformation  energies  in  the  form  of 
(controlled-continuity)  generalized  splines  [75].  The  controlled-continuity  spline  is  a 
generalization  of  Tikhonov  and  Arsenin  stabilizer  [80],  and  can  formally  be  regarded 
as  regularizing  the  ill-posed  problems  [56],  by  recasting  them  as  well-posed  functional 
minimization  problems. 

Active  or  deformable  model  can  represent  broad  classes  of  shapes  by  employing 
geometric  representations  that  involve  many  degrees  of  freedom,  such  as  splines. 
The  degrees  of  freedom  are  generally  not  permitted  to  evolve  independently,  but  are 
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governed  by  physical  principles  that  bestow  intuitively  meaningful  behavior.  The 
use  of  physical  principles,  especially  elasticity  theory,  generally  within  a  Lagrangian 
setting,  give  rise  to  the  name  "deformable  models"  or  "active  models."  The  physical 
interpretation  views  deformable  models  as  elastic  bodies  which  respond  naturally  to 
applied  forces  and  constraints.  Typically,  deformation  energy  functions  are  defined  in 
terms  of  the  geometric  degrees  of  freedom  for  the  deformable  models  as  for  physically 
elastic  objects.  The  energy  of  the  model  grows  monotonically  as  it  deforms  away 
from  a  specified  natural  or  "rest"  shape  and  often  includes  terms  that  constrain 
the  smoothness  or  symmetry  of  the  model.  In  the  Lagrangian  setting,  the  elastic 
forces  internal  to  the  model  contribute  to  the  deformation  energy,  while  the  external 
potential  energy  functions  are  defined  in  terms  of  the  data  of  interest  to  which  the 
model  is  to  be  fitted.  These  potential  energies  give  rise  to  external  forces  which 
deform  the  model  such  that  it  maximally  fits  the  data. 

From  a  theoretical  point  of  view,  the  concept  of  active  or  deformable  models 
reflects  the  integrated  fluency  of  geometry,  physics,  and  approximation  theory.  Geom- 
etry serves  to  represent  object  shape,  physics  imposes  constraints  on  how  the  shape 
may  vary  over  space  and  time,  and  the  theory  of  optimal  approximation  provides  the 
formal  support  to  the  mechanisms  for  fitting  the  models  to  measured  data. 

The  active/deformable  model  that  has  attracted  the  most  attention  in  the  com- 
puter vision  and  computer  graphics  community  is  popularly  known  as  "snakes"  [5,  31]. 
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Snakes  or  "active  contour  models"  represent  a  special  case  of  the  general  multidimen- 
sional deformable  model  theory.  We  will  review  their  formulation  in  this  chapter  in  or- 
der to  illustrate  the  basic  mathematical  foundation  of  the  controlling  curves/surfaces 
in  our  proposed  shape  modeling  scheme,  namely,  the  snake  pedal  model. 

*  2.1    Snakes:  Energy  Minimizing  Active  Contours 

Snakes  are  energy  minimizing  splines  guided  by  internal  and  external  forces. 
Internal  forces  act  as  smoothness  constraints,  while  external  forces  pull  the  snake 
toward  the  desired  features  such  as  edges  and  image  boundaries.  By  initializing  the 
snake  around  the  target  object,  it  can  find  the  desired  object  contours  by  minimizing  a 
user-defined  energy  function.  Because  the  way  the  contours  slither  while  minimizing 
their  energy,  we  call  them  "snakes."  The  ability  of  the  snakes  to  converge  to  the 
outline  of  the  object  despite  the  intensity  gradient  inside  the  object  makes  them  very 
useful  in  many  visual  tasks,  such  as  image  segmentation,  boundary  detection,  motion 
tracking,  and  stereo  matching. 

Geometrically,  a  snake  consists  of  a  set  of  discrete  points,  effectively  connected 
by  straight  lines.  Each  point  has  a  position,  given  by  the  coordinates  in  the  plane 
{x,y),  and  a  snake  is  completely  specified  by  the  number  and  coordinates  of  these 
points,  namely,  p(s)  =  {x{s),y{s))  with  the  domain  parameter  s  e  [0, 1].  When  the 
snake  is  a  closed  curve,  we  apply  the  periodic  boundary  condition,  p(0)  =  p(l),  on 
the  snake. 

A  snake  is  controlled  by  minimizing  a  functional  which  converts  high-level  con- 
tour information,  such  as  curvature  and  discontinuity,  and  low-level Arsenin  image 
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information,  such  as  edge  gradients  into  an  energy  represented  by 

E{s)  =  r  E{p{s))  ds  =  [\eUp{s))  +  EeMs))}  ds  (2.1) 
Jo  Jo 

where,  Eint  and  Ee^t  are  the  internal  and  external  energies  associated  with  the  snake. 
The  shape  and  position  of  the  snake  corresponds  to  the  minimum  of  this  functional. 
The  first  term  of  the  above  functional 


EinMs))  =  Wi{s) 


dp 

ds 


d'^p 


ds' 


(2.2) 


is  the  internal  deformation  energy  of  the  snake.  It  characterizes  the  deformation 
of  an  elastic  contour,  Wi{s)  and  W2{s)  are  non-negative  parameters  controlling  the 
"tension"  and  "rigidity"  of  this  elastic  contour. 

The  second  term  of  (2.1)  consists  of  the  external  potential  that  guides  the  snake 
toward  desired  features.  One  form  of  this  potential  is  a  potential  function  defined 
on  the  image  plane,  whose  local  minima  coincide  with  intensity  extrema,  edges,  and 
other  image  features  of  interest.  For  example,  if  we  want  to  attract  the  contour 
toward  the  edges  of  a  gray  level  image  I{x,  y),  the  potential  takes  the  form 


EexMs))  =        Po{p{s))  ds, 

Jo 


(2.3) 


with 


Po{p{s))  =  Po{x,  y)  =  -c  I VG,  *  I{x,  y)\ , 


(2.4) 
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where  c  controls  the  magnitude  of  the  potential,  V  is  the  gradient  operator,  and 
Ga  *  I  denotes  the  image  convolved  with  a  Gaussian  smoothing  filter,  and  a  is  the 
width  of  the  filter  which  controls  the  spatial  extent  of  the  local  minima  of  Eext- 

Another  form  of  the  external  energy  potential  is  the  user  defined  external  con- 
straint potential  such  as  interactive  springs,  anchored  springs,  and  volcanos  [31]. 
For  example,  the  snake  can  be  pulled  in  the  direction  of  the  mouse  cursor  location 
i^m,  Urn)  by  usiug  a  spring  potential  Ee^t  =  c[{x  —  +  {y  —  J/m)^],  where  c  controls 
the  stiffness  of  the  spring.  Note  that  the  points  on  the  snake  that  are  affected  by  the 
spring  can  be  restricted  to  a  small  section  of  the  snake  closest  to  the  mouse. 

The  combination  of  external  image  potentials  and  external  constraint  potentials 
allows  snakes  to  extract  and  represent  a  broad  spectrum  of  shapes.  Furthermore, 
external  constraint  potentials  are  an  effective,  flexible  means  from  which  high-level 
control  mechanisms  can  guide  shape  recovery,  forming  a  sound  basis  for  fully  auto- 
matic image  analysis. 

Using  the  calculus  of  variations,  we  obtain  the  contours  p(s)  that  minimize  the 
energy  E.  They  must  satisfy  the  Euler-Lagrange  equation 

This  vector-valued  partial  differential  equation  expresses  the  balance  of  internal  and 
external  forces  when  the  contour  rests  at  equilibrium.  The  first  two  terms  represent 
the  internal  stretching  and  bending  forces  respectively,  while  the  third  term  represents 
the  external  forces  that  guide  the  snake  to  the  image  data. 
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In  order  to  compute  numerically  a  minimum  energy  solution,  it  is  necessary  to 
discretize  the  energy  (2.1).  By  using  the  finite  element  method  [87]  or  the  finite 
difference  method  [70],  the  discrete  form  of  the  energy  E{p{s))  can  be  written  as: 


where  K  is  the  stiffness  matrix,  and  E{p)  is  the  discrete  version  of  the  external 
potential.  Differentiating  (2.6)  with  respect  to  p  results  in  the  minimum  energy 
solution,  which  is  equivalent  to  solving  the  set  of  nonlinear  algebraic  equations: 


where  f{x,y)  =  {dPo/dx,dPo/dy)  is  the  generalized  external  force  vector. 

It  is  possible  to  introduce  dynamics  into  the  above  static  formulation  of  the 
snake.  This  may  be  done  by  introducing  a  dissipation  energy  functional  which  dis- 
sipates the  kinetic  energy  during  the  snakes  motion.  The  motion  equation  poses  an 
initial  boundary  value  problem  which  can  be  solved  using  efficient  numerical  tech- 
niques leading  to  a  real  time  response  from  the  snake,  to  applied  forces.  The  initial 
value  problem  can  be  solved  using  explicit  Euler  integration.  Let  7  denote  the  Euler 


(time)  step  size,  then  the  equations  for  time  integration  may  be  written  as 


£;(p)  =  ^p^Kp  +  VPo(p) 


(2.6) 


Kp  +  f  (x,  y)  =  0, 


(2.7) 


P 


(K  +  7l)-i(7p*-f(x*,y')). 


(2.8) 
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For  compactness  in  representation,  we  can  replace  the  snake  with  a  B-snake  which  is 
a  snake  constructed  from  B-splines  that  have  built-in  smoothness.  In  this  case,  we 
simply  replace  p(s)  by  XI^o  C'iBi(s),  where  are  the  control  points  of  the  chosen  B- 
spline  basis  and  m  are  the  number  of  control  points.  For  more  details  on  B-snakes, 
we  refer  the  reader  to  Marc  et  al.  [47]. 

2.2    Active  Surface  Models 

A  active/deformable  surface  is  represented  by  a  vector- valued  parametric  rep- 
resentation p{u,v)  =  {x{u,v),y{u,v),  z{u,v))^,  where  {u,v)  are  domain  parameters, 
and  the  parametric  domain  is  the  unit  square  [0, 1]. 

The  deformation  energy  associated  with  the  deformable  surface,  which  simulates 
the  thin-plate-membrane  material  surface,  is  given  by 

eup)  =  I^Mfj ^  (^?) + M(^r + 2(^/ + i'^n  * d..  (2.9) 

Where  Eint  can  be  regarded  as  a  controUed-continuity  spline  defined  by  Terzopoulos 
[75].  As  described  before,  and  wi  and  W2  control  the  tension  and  rigidity  of  the  spline 
respectively.  Increasing  wi  tends  to  decrease  the  surface  area  of  the  material,  while 
increasing  W2  tends  to  make  it  smoother  but  less  flexible.  Fig.  2.1  demonstrates  the 
effect  of  the  value  changes  of  Wi/w2  to  the  material. 

Using  variational  principles,  we  can  obtain  the  Euler-Lagrange  equation  for  (2.9) 
and  its  corresponding  stiflfness  matrix  K.  For  the  details  of  (2.9),  we  refer  the  reader 
to  the  Appendix  A  and  Terzopoulos  and  Metaxas  [78]. 
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(a)  (b)  (c) 


Figure  2.1.  Active  surface  being  pulled  by  a  spring  force  showing  the  effect  of  various 
wi  and  W2  settings:  (a)  Wi  =  0.9,1^2  =  0.1,  {h)wi  =  W2  =  0.5,  and  (c)  Wi  = 
0.01, =  0.99. 


CHAPTER  3 

SNAKE  PEDALS:  GEOMETRIC  MODELS  WITH  PHYSICS-BASED  CONTROL 


In  this  chapter  we  will  first  present  definitions  for  the  regular  pedal  curves/surfaces 
and  illustrate  the  ideas  via  synthesized  examples.  We  will  then  present  our  proposed 
modeling  scheme — "snake  pedal"  model.  Snake  pedals  are  inherently  geometric  mod- 
els, but  they  can  express  local  detail,  as  well  as  cope  with  topology  changes  of  shapes 
automatically.  We  will  illustrate  the  power  of  this  modeling  scheme  via  some  syn- 
thesized examples  in  this  chapter  as  well  as  some  model  fitting  examples  in  Chapter 
8. 

3.1    Definition  of  Regular  Pedal  Curves 

As  shown  in  Fig.  3.1,  let  a  be  a  planar  curve,  the  pedal  curve  of  a  is  defined 
as  the  locus  of  points  on  the  foot  f  of  the  perpendicular  from  a  fixed  point  p  called 
the  control  point  to  a  variable  tangent  of  a.  More  specifically,  let  the  planar  curve 


pedal  point  p 


Figure  3.1.  f  is  on  the  pedal  curve  of  a  with  respect  to  the  pedal  point  p. 


20 


21 


a  be  a  parameterized  curve  with  domain  parameter  t,  let  /3  be  the  pedal  curve  of  at 
with  respect  to  the  fixed  pedal  point  p,  and  let  a{t)  =  g,  /3  (t)  =  f  at  a  particular 
value  of  parameter  t,  the  projection  of  a  (t)  —  p  in  the  direction  Ja  (t)  must  be  P 
(t)  —  p,  where  a'  {t)  is  the  tangent  of  plane  curve  a  (i),  J  :  — y  is  a  linear 
map  given  by 

J{x,y)  =  {-x,y)  (3.1) 

Note  that  applying  the  operation  J  on  a  vector  (x,  y)  in  5?^  Euclidean  space  results 
in  a  new  vector  (— y,  x),  and  J  can  be  geometrically  interpreted  as  a  rotation  by  n/2 
in  a  counterclockwise  direction.  We  can  thus  formally  define  a  regular  pedal  curve 
as  follows  [23]:    ;  ;  .'  ■  . 

Definition  :  The  pedal  curve  of  a  regular  curve  a  :  (c,  d)  — >  with  respect 
to  a  fixed  ( control)  point  p  G  9?^  is  given  by 

,      .  (a(t)  -  p)  ■  J  a ' (t)  ^    ,,.  ,  ^ 

pedal(i)  =  p  +       II  Ja  '(t)  ^^'^^ 

Where  (c,  d)  is  the  parameter  domain  for  t.  The  definitions  of  regular  pedal 
curves  and  surfaces  presented  here  may  be  found  in  any  book  on  elementary  differ- 
ential geometry  text  [23,  38,  42].  Here  we  adopt  the  definition  by  Gray  [23]. 

Some  examples  depicting  the  pedal  curves  of  an  ellipse  for  different  positions 
of  the  pedal  point  are  shown  in  Fig.  3.2.  Usually  the  global  or  core  shape  depicts 
an  overall  size  and  orientation  of  an  object,  and  local  deformations  depict  fine  detail 
in  the  local  positions  of  an  object.  It  is  clear  from  Fig.  3.2  that  the  pedal  curve  is 
capable  of  exhibiting  both  global  and  local  deformations.  For  example,  the  global 
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(a)  (b)  (c)  (d) 

Figure  3.2.  Examples  of  pedal  curves  of  an  ellipse  for  different  pedal  point  positions. 
Pedal  points  are  shown  by  a  dot  in  each  case. 

shape  of  a  pedal  curve  of  an  ellipse  is  approximated  by  a  new  ellipse,  whose  long 
axis  is  almost  the  same  as  that  of  the  original  ellipse,  while  the  short  axis  is  much 
longer  than  that  of  the  original  ellipse.  The  new  ellipse  is  pinched  in  the  direction  of 
the  short  axis  to  produce  the  pedal  curve.  The  location  and  extent  of  this  pinching 
operation  are  determined  by  the  position  of  the  pedal  point.  As  shown  in  Fig.  3.2, 
if  the  pedal  point  is  located  in  the  center  of  the  original  ellipse  (Fig.  3.2  (a)),  the 
pinching  occurs  around  the  center  part  of  the  pedal  curve.  As  the  pedal  point  moves 
to  an  upper  position  (Fig.  3.2  (b)),  the  pinching  position  in  the  pedal  curve  moves 
up  accordingly.  Similarly,  as  the  pedal  point  moves  to  the  right  of  the  center  of  the 
ellipse  (Fig.  3.2  (c)),  more  pinching  occurs  on  the  right  side  in  comparison  to  the  left 
of  the  origin.  This  pinching  operation  constitutes  the  local  deformation  of  the  pedal 
curve,  but  the  type  of  the  local  deformation  that  can  be  exhibited  by  a  pedal  curve  is 
by  no  means  limited  to  the  pinching-type  or  cusp-type  deformations.  More  complex 
local  deformations  can  be  realized  in  the  pedal  curve  as  shown  later  subsequently. 
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In  summary,  by  moving  the  position  of  the  pedal  point,  it  is  possible  to  synthesize 
a  variety  of  global  and  local  deformations  in  regular  pedal  curves/surfaces.  The 
original  curve  a(f)  will  henceforth  be  referred  to  as  the  generator  for  the  pedal 
curve  j3  (t)  and  the  process  of  generating  a  pedal  curve  will  be  referred  to  as  the 
pedaling  operation,  as  depicted  in  Eq.  (3.2).  Note  that  the  word  "generator"  here 
represents  a  curve  in  2D  or  surface  in  3D,  which  differs  from  its  ordinary  meaning  in 
electricity  industry. 

■  3.2    Snake  Pedals:  The  Proposed  Model 

From  the  discussion  in  Section  3.1,  we  observed  that  when  the  pedal  curves 
are  generated  from  ellipses  with  respect  to  single  fixed  pedal  points,  the  shapes  and 
deformations  exhibited  by  regular  pedal  curves  are  quite  limited.  To  overcome  this 
limitation,  we  can  let  the  pedal  points  be  different  for  each  point  of  the  generator 
when  applying  the  pedaling  operation  to  the  generator,  and  thereby  synthesize  more 
general  shapes  of  pedal  curves/surfaces.  We  can  let  these  pedal  points  be  specified 
by  another  curve  p{t),  which  can  be  represented  by  a  active  or  deformable  model 
(described  in  Chapter  2),  that  is,  a  standard  snake  or  a  B-snake  [31,  47]  and  then 
apply  the  pedaling  operation  to  each  point  on  the  generator  a{ti)  with  respect  to  the 
corresponding  pedal  point  p{ti).  Here  we  use  p{t)  to  represent  a  parametric  curve, 
p{ti)  represents  a  particular  point  on  the  curve  p(^)  when  t  =  U.  Note  that  in  this 
generalization  of  regular  pedal  curve  definition,  for  different  points  on  the  generator, 
we  apply  the  pedaling  operation  with  respect  to  different  control  points,  unlike  in  the 
regular  pedal  curves  where  only  a  single  fixed  pedal  point  is  used. 
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snake  pedal 


snake 


(a)  (b) 

Figure  3.3.  (a)  The  process  of  generating  a  snake  pedal  with  an  ellipse  as  a  generator, 
(b)  "snake  pedal"  controlled  by  the  snake  using  an  ellipse  generator. 

In  this  generalization  process,  we  can  employ  a  proper  association  mechanism  to 
establish  a  correspondence  between  the  points  on  the  generator  a(^)  and  the  points 
on  the  curve  p(^).  This  association  scheme  should  keep  the  continuity  property  of  the 
curve  p(t)  and  be  easy  to  implement.  In  accordance  with  these  two  requirements,  we 
developed  a  simple  correspondence  association  mechanism  as  shown  in  Fig.  3.4.  In 
this  mechanism,  for  each  point  (x{ti)  on  the  generator  ,  we  draw  a  straight  line  from 
the  center  of  the  generator  O  to  the  point  a{ti),  and  denote  by  p{ti)  the  intersection 
of  the  line  Oa{ti)  and  the  curve  p{t).  When  applying  the  pedaling  operation  to 
point  a{ti),  we  carry  out  this  operation  with  respect  to  point  p{ti).  The  pedaling 
operation  is  performed  for  each  of  the  points  a{tj),j  =  1, 2, M,  on  the  generator 
resulting  in  a  curve  that  we  call  the  "snake  pedal."  Since  Oa{ti)  and  OpiU)  are 
on  the  same  radial  line  of  the  generator,  we  refer  to  this  association  scheme  as  the 
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0C(t4) 


a(ti)/ 


a(tj) 


"\ a(t2) 


a(ti) 


snake 


generator 


(a) 


Figure  3.4.  Radial- based  correspondence  association  scheme:  each  point  oi{ti)  on  the 
generator  is  associated  with  a  point     on  the  snake. 

radial-based  association  technique.  It  is  obvious  in  this  scheme  that  we  use  the 
same  parameterization  for  the  generator  a(t)  and  curve  p(t)  with  the  same  domain 
parameter  t.  It  can  be  proved  that  snake  pedal  has  the  same  continuity  property  as 
that  of  the  curve  p(t),  given  that  the  generator  is  differentiable,  although  this  may 
not  always  be  the  case. 

The  generator  can  be  either  a  parameterized  or  an  implicit  function  representing 
a  curve.  In  this  chapter,  we  will  focus  on  parameterized  curves  and  in  most  applica- 
tions, we  use  an  ellipse  as  the  generator  in  2D.  If  the  generator  is  an  ellipse,  we  can 
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represent  it  in  a  parametric  form  by 


cos  6  sin  0 
—sin  9  cos  9 


ai  cos  t 
a2  sin  t 


mi 
m2 


(3.3) 


where  ai,a2  are  aspect  ratio  parameters  and  we  use  =  (01,02)^  to  denote  these 
shape  parameters,  9  is  the  rotation  angle  between  the  world  coordinates  and  the 
object  centered  coordinates  which  is  denoted  by  ge  —  (9)'^,  and  gc  =  (mi,  7712)^  is 
the  centroid  of  the  generator  in  the  world  coordinates.  We  collect  all  the  generator 
parameters  into  a  vector  g  =  (gf ,  g^",  gj)"^  =  (ci,  a2,9,  mi,  7712)^,  and  the  vector  g  is 
referred  to  as  the  global  shape  parameter  vector. 

For  the  ellipse,  the  term  Ja'{t)  in  the  Eq.  (3.2)  of  the  pedal  curve  can  be 
expressed  as 


Ja'{t)  = 


ni 

n2 

-ai  sin  t  sin  9  —  a2  cos  t  cos  9 
-ai  sin  t  cos  9  +  a2  cos  t  sin  9 


(3.4) 


Note  that  ||  Ja'lp  =  n^+nl  =  sin'^t+al  cos'^t,  which  is  independent  of  the  rotation 
angle  9.  We  will  see  that  ||Ja:'|p  is  a  very  important  quantity  in  the  model  fitting 
presented  in  later  Chapters. 

The  curve  p{t)  on  which  the  pedal  points  are  located  can  be  of  any  kind,  which 
is  not  limited  to  be  to  a  parameterized  function.  By  using  the  correspondence  asso- 
ciation scheme  described  earlier,  we  can  sample  the  curve  at  the  same  rate  as  that 
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of  the  generator  so  that  every  point  on  the  generator  a  [t)  can  find  a  unique  cor- 
responding point  on  the  snake  p{t).  The  pedaUng  operation  can  then  be  applied  to 
every  point  on  the  generator  with  respect  to  its  corresponding  point  on  the  snake. 
We  use  a  snake/B-snake  as  this  controlling  curve  and  get  a  sequence  of  position  vec- 
tors, Pi  =  {xi,yiY ,  by  sampUng  the  snake  curve.  We  collect  all  the  snake  position 
vectors  in  another  vector,  p  =  (pj",  pj, Pm)*^;  where  M  is  the  number  of  sampling 
points.  Vector  p  is  referred  to  as  the  local  shape  parameter  vector.  Vector  p  is 
the  discretized  version  of  continuous  curve  p(t). 

In  our  implementation,  since  the  pedal  points  are  located  on  a  snake  or  B-snake 
and  the  generalized  pedal  curve  is  created  by  applying  the  pedaling  operation  to  the 
generator  with  respect  to  the  points  on  the  snake,  we  dub  this  generalized  pedal 
curve  as  a  snake  pedal.  With  the  radial-based  correspondence  associating  scheme, 
a  snake  pedal  is  completely  determined  by  the  global  shape  parameter,  namely,  the 
generator  parameter  g,  and  the  local  shape  parameter,  namely,  the  snake  position 
vector  p.  Fig.  3.3  demonstrates  the  process  of  creating  the  snake  pedal  curve. 

With  the  correspondence  association  scheme,  the  generalized  pedaling  operation 
can  create  much  broader  classes  of  shapes  with  more  deformations.  We  show  some 
examples  of  snake  pedal  curves  generated  using  snakes  and  an  ellipse  as  the  generator 
in  Fig.  3.5.  Note  the  variety  of  local  deformations  that  can  be  generated  using  this 
modeling  technique.  We  remind  the  reader  that  the  snake  pedal  itself  is  a  geometric 
model  and  that  it  is  not  directly  responsive  to  the  application  of  external  forces  unlike 
the  standard  snake  models  described  in  Chapter  2.   Also,  it  is  easy  to  synthesize 
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(a)  (b)  (c) 


.■      ,  (d)  (e)  (f) 

Figure  3.5.  Examples  of  "snake  pedals"  using  an  ellipse  generator  and  a  snake. 

snake  pedals  whose  topology  is  different  from  that  of  the  generator,  as  will  be  seen  in 
examples  depicted  subsequently.  The  uniqueness  of  our  modeling  scheme  lies  in  the 
ability  to  generate  snake  like  local  deformations  in  a  non-reactive  geometric  shape, 
namely,  the  pedal  curve/surface. 

3.3    Modified  Snake  Pedal  Model 

The  pedal  curve  definition  can  be  modified  slightly  by  subtracting  the  second 
term  from  the  first  term  in  equation  3.2, 

pedal(i)  =  p  -  ^-AJ^>_±lj^  (3.5) 
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(d)  (e)  (f) 

Figure  3.6.  Examples  of  the  modified  "snake  pedals"  using  an  ellipse  generator  and 
a  snake. 

This  modification  allows  us  to  synthesize  a  larger  class  of  shapes  which  are  quite 
different  from  the  shapes  produced  using  Eq.  (3.2).  The  key  feature  of  using  this 
modification  is  that  the  snake  pedal  curve  allows  for  more  local  deformations  including 
shrinkage  and  expansion  of  the  snake  pedal.  The  shrinkage  and  expansion  behavior 
is  controlled  by  the  location  of  the  snake.  When  the  snake  is  located  inside  the 
generator,  the  snake  pedal  exhibits  shrinkage  behavior,  while  expansion  behavior 
appears  when  the  snake  is  outside  the  generator. 

Fig.  3.6  depicts  examples  of  shapes  generated  using  this  modification.  Note  the 
amount  of  local  deformations  allowed  in  the  snake  pedal  curves  using  this  modifica- 
tion. In  Fig.  3.7,  we  compare  the  original  and  "modified"  snake  pedals  using  the 
same  generators,  same  snakes  and  same  correspondence  associating  schemes.  It  is 
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Figure  3.7.  Comparisons  of  original  (red)  and  modified  (blue)  "snake  pedals"  using 
the  same  generators  and  same  snakes  in  both  case. 

obvious  that  the  "modified"  snake  pedals  exhibit  larger  local  deformations  in  both 
Fig.3.7  (a)  and  (b). 

3.4    Handling  Changes  in  Topology 

,  Topology  changes  may  be  achieved  by  applying  the  pedaling  operation  to  a 
single  generator  and  several  pedal  points  (or  snakes).  Fig.  3.8  depicts  examples  of 
the  topology  change  realized  in  2D.  When  the  snake  is  located  inside  the  generator, 
a  single  closed  snake  pedal  curve  is  created  (Fig.  3.8  (a)),  when  the  snake  is  split 
into  two  parts  with  both  parts  located  outside  the  generator,  two  separated  closed 
snake  pedal  curves  are  obtained  (Fig.3.8  (b)).  In  the  second  case,  only  part  of  the 
generator  is  needed  to  produce  the  two  snake  pedal  curves.  Further  study  shows  that 
any  number  of  closed  snake  pedals  can  be  generated  from  one  generator  by  using  this 
procedure,  if  the  position  of  snakes  and  parts  of  the  generator  are  properly  chosen. 
When  using  the  correspondence  associating  scheme  discussed  previously,  the  number 
of  closed  snake  pedal  contours  are  determined  by  the  number  of  closed  snakes.  More 
generally,  the  topological  property  of  the  snake  pedal  is  determined  by  the  topology  of 
the  controlling  snake.  Therefore,  we  can  accomplish  the  automatic  topological  change 
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(a)  (b)  (c)  (d) 

Figure  3.8.  Examples  of  topology  change,  (a)  depicts  a  single  snake  pedal  curve 
(solid)  obtained  from  an  ellipse  generator  (dashed)  and  a  single  pedal  point,  (b) 
same  as  (a)  except  two  distinct  pedal  points  located  outside  the  generator  are  used, 
(c)  same  as  (a)  except  several  pedal  points  on  a  snake  are  used,  (d)  same  as  (c) 
except  a  discontinuous  snake  is  used. 

for  the  snake  pedal  model  if  the  topology  of  the  snake  can  be  updated  automatically 
in  the  data  fitting  applications.  In  Chapter  6,  we  will  introduce  a  novel  level-set  based 
approach  to  realize  the  automatic  topological  change  for  the  snake  pedal  model. 

3.5    Pedal  Surfaces 

A  regular  pedal  surface  is  the  surface  analog  of  the  pedal  curve.  It  is  the  locus 
of  the  points  on  the  foot  of  the  perpendicular  from  a  fixed  pedal  point  to  a  variable 
tangent  plane  of  the  surface.  More  precisely,  it  is  defined  as  follows  [23]. 

Definition  :  The  pedal  surface  of  a  parameterized  surface  or  a  patch  ct  : 
U  — >  3?^  with  respect  to  a  point  p  e      is  defined  by 
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Where  U  is  an  open  subset  of  3?^,  cx{u,v)  is  a  parameterized  generator  surface 
with  parameters  u  and  v  and  ttu  and  a„  are  tangent  vectors  in  the  u,  v  directions. 
For  an  ellipsoid  generator,  we  have 


Qi  cos  U  COS  V 

mi 

a{u,  v)  =  H 

02  COS  u  sin  V 

+ 

m2 

aa  sin  v 

(3.7) 


Where  ai ,  02  and  03  are  aspect  ratio  parameters  which  can  be  collected  in  gs  = 
(oi,  02,  as)^.  We  use  gc  =  (mi,  7712,  ma) ^  denote  the  centroid  of  the  generator  in  the 
world  coordinates.  R  is  the  rotation  matrix  between  the  world  coordinates  and  the 
object  centered  coordinates.  We  use  quaternions  [63]  to  represent  the  rotation  matrix 
because  the  rotation  matrix  is  easier  to  update  through  the  quaternion  update  and 
remains  orthogonal  in  the  mean  time.  A  quaternion  [s,  v]  defines  a  rotation  of  the 
model  from  its  reference  position  through  an  angle  9  =  2  cos~^s  around  an  axis 
aligned  with  vector  v  =  {vi,V2,V3)'^ .  The  rotation  matrix  corresponds  to  [s,v]  is 

1  -  2{vl  +  vj)     2{viV2  -  SV3)       2{ViV3  +  SV2) 

=    2{viV2  +  SV3)     1  -  2{vj  +  vj)   2{v2V3  -  svi)      ■  (3-8) 

2{ViV3  -  SV2       2{v2Vs  +  SVi)      1  -  2{vj  +  vl) 

We  collect  the  quaternion  parameters  into  a  new  vector  g^  =  (s,  q^)  =  (s,  Vi,V2,V3) 
The  ellipsoid  generator  is  then  completely  represented  by  g  =  (gj',  gj,  gj)^.  As  in 
2D  case,  g  is  referred  to  as  the  global  shape  parameter  vector. 
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(d)  (e)  (f) 

't       '        N  1 

Figure  3.9.  Examples  of  pedal  surfaces  with  fixed  pedal  points  located  inside  ellip- 
soidal generator  surfaces. 

Fig.  3.9  depicts  examples  of  pedal  surfaces  synthesized  using  an  ellipsoidal 
generator  surface  and  a  fixed  pedal  point  located  inside  the  generator  surface  in  each 
example.  Note  that  change  in  topology  is  achieved  by  simply  choosing  an  appropriate 
location  of  the  pedal  point  as  shown  in  Fig.  3.9.  However,  this  type  of  topology 
change  is  diflferent  from  what  we  described  in  Fig.  3.8. 

As  in  the  2D  case,  we  can  let  the  pedal  point  vary  for  each  point  on  the  generator 
surface.  Thus  we  have  the  snake  pedal  surfaces  in  3D  whose  shape  can  be  controlled 
by  snakes  which  are  either  curves  or  surfaces  in  3D.  Fig.  3.10  depicts  some  examples  of 
snake  pedal  surfaces  in  3D.  They  are  synthesized  by  using  an  ellipsoid  as  a  generator 
and  a  snake  curve  in  3D  which  lies  inside  the  generator.  The  curve  is  shown  outside 
of  the  generator  adjacent  to  the  snake  pedal  model  for  the  clarity  of  visualization. 
Note  that  the  snake  pedal  shape  can  be  manipulated  by  applying  external  forces  to 
the  snake  as  illustrated  in  the  Fig.  3.10  (a)-(d). 
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(c)  (d) 

Figure  3.10.  Examples  of  shaping  snake  pedal  surfaces  with  snakes.  Each  figure  was 
generated  by  applying  interactive  mouse-based  forces  to  the  snake. 

For  the  3D  snake  pedal  examples  obtained  by  ellipsoids  and  snake  surfaces,  we 
refer  the  reader  to  the  next  chapters. 

3.5.1    Uniform  Tessellation  of  the  Model 

From  the  definition  (3.6)  of  the  pedal  surface,  the  points  on  the  pedal  surface 
must  be  located  in  the  tangent  of  the  generator  for  every  point  on  the  generator 
specified  by  {u,v).  For  an  ellipsoid  generator,  the  evenly  spaced  domain  parameters 
u  and  V  are  often  referred  to  as  the  "latitude"  and  "longitude."  It  should  be  noted 
that  the  snake  pedal  model  exhibits  nonuniformly  tessellated  surface  in  spite  of  the 
domain  parameters  u  and  v  being  uniformly  tessellated/sampled.  This  is  obvious 
from  the  nonlinear  function  defining  the  pedal  surface  in  Eq.  (3.6).  For  example,  the 
tangent  of  the  generator  undergoes  a  rapid  change  (in  orientation)  as  u      n/2,  and 
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a  slow  change  (in  orientation)  as  u  0.  Therefore,  for  an  evenly  sampled  generator 
and  snake,  when  we  use  the  association  correspondence  mechanism  to  produce  a 
snake  pedal  surface,  the  snake  pedal  model  will  not  be  uniformly  tessellated. 

This  uneven  tessellation  phenomenon  is  not  noticeable  when  the  aspect  ratio 
parameters  Oi ,  02  and  03  of  the  ellipsoid  generator  do  not  differ  much  from  each 
other.  But  when  the  difference  of  the  aspect  ratio  parameters  increases,  this  non- 
uniformity  becomes  obvious,  and  leads  to  the  distortion  of  the  deformation  properties 
of  the  snake  pedal  surfaces.  In  addition,  numerical  instabilities  may  be  caused  in 
the  computation  of  derivatives  for  nonuniformly  tessellated  snake  pedal.  Hence  a 
uniform  tessellation  scheme  for  pedal  surfaces  must  be  employed  to  overcome  this 
problem.  It  is  a  natural  way  to  first  construct  uniformly  spaced  tangent  directions  in 
the  parameter  direction  u,  then  apply  the  same  technique  to  achieve  uniform  spacing 
along  the  parameter  direction  v.  We  now  present  our  scheme  for  the  construction  of 
uniform  tessellation  of  snake  pedal  surfaces. 

When  constructing  the  uniform  tessellation  of  an  ellipsoid  generator  in  the  pa- 
rameter direction  u,  we  actually  produce  a  uniformity  in  the  tangent  direction  for 
an  ellipse  along  a  fixed  longitude  {v  is  fixed  in  this  case).  Since  we  do  not  need  to 
consider  the  rotation  and  translation  of  the  ellipse  when  constructing  the  uniform 
tessellation,  for  an  arbitrary  point  (xo,yo)  on  an  ellipse,  the  standard  parametric 
representation  is: 


V. 


(3.9) 
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The  straight  line  equation  for  the  tangent  of  this  elUpse  on  the  point  (xq,  yo)  is 


^  +  ^  =  1,  (3.10) 


Substituting  Eq.  (3.9)  into  (3.10),  we  obtain 


a2  COSU  a2  in  ^^\ 

y  =  : — x  +  - — .  (3.11) 

ai  sinu  smu 


The  slope  of  this  tangent  line  is:  —  ^ ,  which  corresponds  to  an  angle  of  0  = 

tan~^{  ^  cotu).  The  angle  9  is  not  uniformly  spaced  when  the  variable  u  is  evenly 
spaced,  which  results  in  a  non-uniform  tessellation  of  the  model.  The  objective  of 
the  uniform  tessellation  is  to  modify  the  standard  ellipse  parametric  formula  so  as  to 
construct  a  uniformly  tessellated  9  when  u  is  evenly  spaced.  We  change  the  standard 
ellipse  equation  (3.9)  to  : 

xq   =   ai  cos  f{u) 

(3.12) 

yo    =  a2sinf{u), 

where  f{u)  is  the  function  to  be  solved.  The  angle  of  the  tangent  direction  of  the 
ellipse  accordingly  becomes  9  =  tan~^{  ^  cot  f{u) ).  To  make  this  angle  uniformly 
tessellated,  we  must  set  .  ^ 


^{tan-\-"^cotf{u))}  =  C,  (3.13) 


37 


where  C  is  a  constant.  The  boundary  condition  is 


fin) 


2  ' 


(3.14) 


< 


u  =  + 


2 


/(«) 


Basic  algebraic  manipulation  of  Eq.  (3.13)  results  in  the  following  ordinary  differen- 
tial equation  (ODE)  for  f{u) 


where  f'{u)  is  the  derivative  of  f{u)  with  respect  to  u.  This  ordinary  differential 
equation  can  be  solved  numerically.  We  use  a  third  order  Runge-Kutta  method.  The 
result  is  shown  in  Fig.  3.11  (a).  It  can  be  seen  that  when  u  is  small  {u  — >  0), 
f{u)  changes  much  faster  than  u,  which  tends  to  make  up  the  fact  that  the  tangent 
direction  of  the  ellipse  changes  slowly  in  this  region.  While  when  u  — 7r/2,  f{u) 
changes  much  slower  than  u,  which  is  in  an  attempt  to  slow  down  the  rapid  change 
in  tangent  direction  in  this  area.  After  changing  u  to  f{u)  in  the  ellipse  expression, 
the  change  in  tangent  directions  will  become  even.  Carrying  out  the  same  procedure 
on  the  parameter  v,  we  obtain  the  f{v)  ~  v  relationship  as  shown  in  Fig.  3.11  (b). 
Note  that  the  retessellation  scheme  needs  to  be  done  whenever  the  parameters  of  the 
snake  pedal  are  subject  to  large  changes. 


f'{u)  =  C  ■  {'^szn'f{u)  +  ^cos'f{u)), 


(3.15) 
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(a)  (b) 

'  "  r 

Figure  3.11.  Uniform  tessellation  result:  (a)  f{u)  ~  u  relation,  (b)  f{v)  ~  v  relation. 

3.6    Summary:  Salient  Features  of  the  Model 

Compared  with  the  hybrid  models  mention  in  Section  1.2.1,  the  snake  pedal 
models  are  bestowed  with  several  unique  features  which  make  them  very  powerful 
and  attractive  for  use  in  shape  synthesis  and  recovery  applications.  A  list  of  these 
features  are: 

1.  The  snake  pedal  model  is  inherently  geometric  but  can  exhibit  both  local  and 
global  deformations  as  would  a  physics-based  model. 

2.  The  model  can  exhibit  large  global  deformations  such  as  bending,  twisting,  and 
tapering  without  explicitly  incorporating  parameters  for  the  same. 

3.  The  modeling  scheme  allows  (for  the  first  time)  the  concept  of  a  global  shape 
to  be  incorporated  into  the  curve/surface  evolution  framework. 


4.  The  model  permits  a  level-set  implementation  and  thereby  seamlessly  allows 
for  topological  changes  in  shape  during  evolution. 


CHAPTER  4 

PROBABILISTIC  FRAMEWORK  AND  MODEL  LEARNING 


In  Chapter  3,  we  described  our  snake  pedal  model,  which  is  fully  characterized 
by  the  shape  parameter  vector  q  =  (g^,p^)^  with  g,p  being  the  global  and  local 
shape  parameter  vectors  respectively.  In  many  applications,  the  general  shape,  loca- 
tion and  orientation  of  objects  is  often  known  a  priori  and  this  knowledge  may  be 
incorporated  in  the  form  of  initial  conditions,  data  constraints,  model  shape  param- 
eter constraints,  or  into  the  model  fitting  procedures.  For  example,  in  most  medical 
images,  the  boundaries  of  the  anatomical  structures  of  interest  are  seldom  well  de- 
lineated from  the  background,  either  because  the  image  quality  is  poor  or  the  shape 
themselves  gradually  blend  into  the  surroundings.  In  many  medical  imaging  applica- 
tions, it  is  necessary  to  use  prior  knowledge  about  the  shape  of  interest  by  learning 
the  shape  from  a  human  expert  or  from  a  shape  atlas.  Proper  learning  schemes  must 
be  developed  to  incorporate  the  human  expert  knowledge  or  the  knowledge  from  an 
atlas.  In  addition,  for  automatic  image  interpretation,  it  is  necessary  to  have  a  model 
that  not  only  describes  the  size,  shape,  location  and  orientation  of  the  target  shape 
but  that  also  permits  expected  variations  in  these  characteristics. 

To  meet  above  requirements,  a  number  of  researchers  have  cast  the  model  fitting 
process  in  a  probabilistic  framework  to  include  prior  knowledge  of  object  shape  by 
incorporating  prior  probability  distributions  on  the  shape  variables  to  be  estimated 
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[29,  69,  81,  86].  In  Staib  and  Duncan  [69],  a  deformable  contour  model  was  used 
on  2D  echocardiograms  and  MR  images  to  extract  the  left  ventricle  (LV)  of  the 
heart  and  the  corpus  callosum  of  the  brain  respectively.  This  closed  contour  model  is 
parameterized  using  an  elliptic  Fourier  decomposition  and  a  priori  shape  information 
is  included  as  a  spatial  probability  expressed  through  the  likelihood  of  each  model 
parameter.  The  model  parameter  probability  distributions  are  derived  from  a  set  of 
example  object  boundaries  and  serve  to  bias  the  contour  model  towards  expected  or 
more  likely  shapes. 

Szekely  et  al.  [71]  have  also  developed  Fourier  parameterized  models.  In  ad- 
dition, they  have  added  elasticity  to  their  models  to  create  "Fourier  snakes"  in  2D 
and  elastically  deformable  Fourier  surface  models  in  3D.  By  using  the  Fourier  pa- 
rameterization followed  by  a  statistical  analysis  of  a  training  set,  they  define  mean 
organ  models  and  their  eigen-deformations.  An  elastic  fit  of  the  mean  model  in  the 
subspace  of  eigen-modes  restricts  possible  deformations  and  finds  an  optimal  match 
between  the  model  surface  and  boundary  candidates. 

Cootes  et  al.  [13]  and  Hill  et  al.  [28]  present  a  statistically  based  technique 
for  building  deformable  shape  templates  and  use  these  models  to  segment  various 
anatomical  structures  from  2D  and  3D  medical  images.  The  statistical  parameteriza- 
tion provides  global  shape  constraints  and  allows  the  model  to  deform  only  in  ways 
implied  by  the  training  set.  These  shape  models  represent  anatomical  structures  by 
sets  of  landmark  points  which  are  placed  in  the  same  way  on  an  object  boundary 
in  each  input  image.  For  example,  to  extract  the  LV  from  echocardiograms,  they 
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choose  points  around  the  ventricle  boundary,  the  nearby  edge  of  the  right  ventricle, 
and  the  top  of  the  left  atrium.  These  points  can  be  connected  to  form  a  deformable 
contour.  By  examining  the  statistics  of  training  sets  of  hand-labeled  medical  images, 
and  using  principal  component  analysis,  a  shape  model  is  derived  that  described  the 
average  positions  and  the  major  modes  of  variation  of  the  points.  New  shapes  are 
generated  using  the  mean  shape  and  a  weighted  sum  of  the  major  modes  of  varia- 
tion. Shape  boundaries  are  then  segmented  using  this  "point  distribution  model"  by 
examining  a  region  around  each  model  point  to  calculate  the  displacement  required 
to  move  it  towards  the  boundary,  these  displacements  are  then  used  to  update  the 
shape  parameter  weights. 

In  this  thesis,  we  cast  the  model  fitting  process  in  a  Bayesian  framework  and 
incorporate  prior  distributions  on  the  vector  of  shape  parameters  that  are  being  es- 
timated. By  using  a  probabilistic  approach  rather  than  a  deterministic  approach, 
we  can  develop  probabilistic  models  of  data  acquisition  which  closely  match  the 
characteristics  of  real  image  sensors,  and  more  importantly,  we  can  determine  the 
uncertainty  in  our  estimate,  which  can  be  used  (if  necessary)  when  integrating  mul- 
tiple sources  of  data/information  [73].  We  now  present  the  details  of  the  Bayesian 

probabilistic  framework. 

f  ■■  -- 

4.1    Prior  Model  in  the  Bavesian  Framework 

As  discussed  earlier,  in  many  applications  such  as  medical  image  analysis,  the 
boundaries  of  anatomical  structures  are  usually  "fuzzy"  because  of  the  partial  volume 
averaging  effects  and  other  effects.  We  must  incorporate  the  prior  information  about 
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the  shape  of  interest  into  the  shape  representation.  Therefore,  in  the  prior  model,  in 
addition  to  generic  smoothness  assumptions,  we  should  also  incorporate  into  the  prior 
model  the  information  about  the  possible  shapes  that  an  observed  object  can  assume. 
We  express  the  probability  of  a  given  configuration  q  through  a  Gibbs  distribution 
[20],  which  assigns  a  low  probability  to  configurations  with  large  energy  and  vise 
versa.  The  configuration  corresponding  to  Ep  =  0  represents  the  rest  state  of  the 
model.  In  the  Gibbs  distribution,  the  internal  deformation  energy,  which  is  similar 
to  the  internal  deformation  energy  in  generalized  splines  described  in  Chapter  2,  is 
transfered  into  a  probabilistic  distribution  by: 

p(q)  =  ^exp{-Ep{q)),  (4.1) 

where  Zp  is  a  normalizing  factor  and  the  energy  Ep{-)  governs  the  deformation  of 
the  model  away  from  its  natural  shape.  There  are  two  parts  contained  in  this  defor- 
mation energy,  one  part  corresponds  to  the  local  shape  deformation,  and  the  other 
part  corresponds  to  the  global  deformation.  For  the  local  shape  deformation,  if  we 
denote  the  rest  shape  of  the  local  representation  (snake  positions)  by  p,  which  can  be 
obtained  by  the  scheme  discussed  in  Section  4.4,  the  internal  energy  in  a  continuous 
form  is  given  by 
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with  the  discrete  version 


l(p_p)K(p_p),  (4.3) 


where  K  is  the  stiffness  matrix  for  the  membrane  plus  thin  plate  deformation  energy 
imposed  on  the  snake.  Note  that  Eq.  (4.2)  and  (2.9)  have  the  similar  form,  if  we  set 
p  =  0,  Eq.  (4.2)  goes  to  (2.9).  Eq.  (2.9)  is  the  special  case  of  (4.2)  where  the  local 
rest  shape  is  zero. 

For  the  global  shape  representation,  if  no  additional  information  is  made  about 
the  shape  being  modeled,  the  underlying  global  shape  of  the  model,  the  generator, 
can  take  the  mean  value  of  the  shapes  that  has  been  learned.  By  extending  the  above 
internal  deformation  energy,  we  can  include  the  global  shape  parameters  of  the  model 
in  our  prior  statistics  in  the  similar  way.  The  extended  energy  can  be  expressed  as: 

^      _     ^{local)  _|_  ^(global) 

=  lip-pr  K  (p-p) 

+   l(gs-gs)^  C;'  (gs-gs)  (4-4) 

+   l(gc-gc)^  (gc-gc). 

Where  gs,  ge,  and  gc  are  the  rest  state  of  the  aspect  ratio  vector,  rotation  vector,  and 
translation  vector  of  the  generator  discussed  in  Chapter  3.  Cs''\C0~^,  and  Cc"^ 
represent  the  covariance  matrices  which  are  accumulated  in  the  iterative  training 
phase  described  in  Section  4.4. 
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4.2    Sensor  Model  in  the  Bayesian  Framework 

The  data  acquisition  (sensor)  model  in  the  Bayesian  framework  is  based  on 
linear  measurements  with  Gaussian  noise  and  given  by: 


p{B/q)  =  ^exp{-ED{D,q)),  (4.5) 


where  Zd  is  a  constant,  D  is  the  data,  and  jE'd(-)  can  be  either  the  edge-based 
potential  energy  synthesized  from  image  data,  or  the  potential  energy  in  the  springs 
attached  from  3D  data  points  to  the  surface  of  the  snake  pedal,  constraining  it  to 
conform  to  the  observed  data.  These  two  types  of  sensor  model  energies  are  expressed 
as: 

{h  Et  c{Di  -  Xi)2        spring  energy, 
(4.6) 
— c||(Gct  *  V/(x))|p   image  potential. 

In  the  first  form  of  energy  Ed{-),  c  represents  the  uncertainty  of  the  sensor  measure- 
ments and  Xj  is  the  closest  model  point  from  the  given  data  Di.  In  the  second  form  of 
energy  -E'd(-))  ^^(x)  is  the  gradient  of  the  image  intensity  /(x),  is  the  Gaussian 
filter  mask  with  width  a,  and  *  represents  the  convolution  operation. 

Taking  the  partial  derivative  of      with  respect  to  q  yields  the  external  force: 
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The  external  force  pulls  the  model  toward  the  image  boundary  or  range  data,  thus 
molding  the  model  into  the  desired  shapes.  This  process  is  known  as  the  model  fitting 
process. 

4.3    Maximum  a  Posterior 

In  the  probabilistic  framework,  combining  the  prior  model  (4.1)  and  sensor 
model  (4.5)  by  using  Bayes  rule,  we  obtain  the  posterior  distribution 

^(^ID)  .  Pj^m^  ^  i.xp(-£(q)),  (4.8) 

where  Z  is  a  constant,  and  the  total  energy  is: 

E{ci)  =  Ep{q)  +  ED{B,q).  (4.9) 

In  the  Bayesian  framework,  the  objective  is  to  compute  the  Maximum  a  Posterior 
(MAP)  estimate,  that  is,  the  value  of  q  that  maximizes  the  conditional  probability 
p(q|D).  This  estimation  provides  the  same  result  as  finding  the  minimum  energy  con- 
figuration of  the  physically  based  model.  However,  compared  with  the  deterministic 
model,  the  probabilistic  framework  has  several  advantages  as  described  earlier. 

4.4    Supervised  Learning  Scheme 

Once  we  have  established  the  Bayesian  statistical  modeling  scheme,  we  propose  a 
supervised  learning  scheme  to  incorporate  the  information  about  the  shape  of  interest 
to  the  model.  In  this  learning  scheme,  we  divide  the  shape  recovery  process  into  a 
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training  phase  and  an  operational  phase.  In  the  training  phase,  prior  knowledge  about 
the  shapes  to  be  recovered  from  the  image  data  is  collected  via  a  supervised  learning 
technique.  A  training  sequence  of  several  iterations  provides  samples  of  correctly 
reconstructed  shapes  for  specific  classes  of  shapes,  for  example,  the  caudate  nucleus 
and  the  hippocampus  from  the  human  brain.  We  use  Pg  denote  the  semi-global  part 
our  model,  which  may  correspond  to  the  coarse  level  wavelet  coefficients  of  the  local 
shape  parameter  vector  p.  The  parameter  sets  q(z)  corresponding  to  these  shapes 
represent  the  training  sequence  for  the  global  part  g  and  semi-global  part  Pg  of  our 
prior  model  as  shown  in  Fig.  (4.1).  With  each  iteration  i,  we  improve  estimates  of  the 
mean  vector  q  and  covariance  matrix  C  through  a  simple  equation  for  accumulating 
measurements  as  indicated  in  Fig.  (4.1).  Since  we  update  only  the  statistics  of  the 
global  parameters  of  the  model,  the  local  deformation  remains  constrained  only  by 
smoothness  assumptions  captured  in  the  stiffness  matrix  K.  Therefore,  most  of  the 
flexibility  of  the  locally  deformable  surface  remains  unchanged. 

In  the  operational  phase,  the  mean  and  covariance  of  the  model  state  q  are 
used  in  initializing  the  model  for  fitting  the  model  to  a  data  set  not  included  in  the 
training  phase.  The  main  difference  between  training  and  operational  phase  is  that, 
in  the  former,  the  model  initialization  for  the  fitting  is  far  from  the  final  desired  fit 
because  no  prior  information  is  available.  Most  existing  schemes  in  literature  that  use 
a  priori  knowledge  about  anatomical  shapes  of  interest  manually  construct  their  prior 
models  [69,  73]  which  can  be  a  very  tedious  process.  The  model  fitting  scheme  that 
was  described  in  this  thesis  in  conjunction  with  the  supervised  learning  mentioned 


Set  of  training  data  samples 
D(z) 

obtained  from  sensor 


Surface  fitting  by  maximizing 
p(q|D(i))    ,  where  current  prior 
model  Pi^)    is  used 

q(0  -  final  result 

Iterative  improvment  of  the  mean 
and  covariance  of  the  global  model 

q{i)  =  iq(i)  +  i^q(z  -  1) 


=  fqWq^W  +  ^m2(7  -  1) 
Cii)  =  m,{z)  -  q(z)q^(i) 


Figure  4.1  Supervised  learning  to  incorporate  information  into  the  prior  model 
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earlier  and  discussed  in  detail  by  Vemuri  and  Radisavljevic  [81]  can  be  used  for  fast 
automatic  construction  of  prior  models.  We  will  devote  the  next  chapter  to  present 
the  fast  numerical  algorithms  for  automatic  model  construction  via  fitting. 


CHAPTER  5 

NUMERICAL  SOLUTIONS  FOR  MODEL  FITTING 
WITH  SNAKE  PEDAL  MODELS 


In  this  chapter  we  will  develop  numerical  techniques  for  fitting  our  snake  pedal 
model  to  a  sparse  set  of  data  points  in  2D/3D  as  well  as  to  image  data  for  shape 
recovery.  The  objective  of  model  fitting  is  to  mold  or  sculpt  the  snake  pedal  model 
so  that  the  model  passes  through  the  set  of  data  points  or  conforms  to  the  de- 
sired image  boundaries  as  closely  as  possible,  while  maintaining  a  certain  degree  of 
smoothness.  This  fitting  process  is  equivalent  to  finding  the  minimum  of  the  to- 
tal energy  (4.9)  of  the  model  described  in  Chapter  4,  which  consists  of  two  terms, 
namely,  the  internal  strain  energy  Ep{-)  corresponding  to  the  snake  deformation  en- 
ergy and  the  external  energy  Ed{-).  The  energy  minimization  is  posed  as  a  non- 
linear optimization  problem  and  can  be  solved  numerically.  There  are  two  types  of 
parameters  which  need  to  be  determined  in  this  optimization  process.  One  is  the 
global  shape  parameters,  namely,  generator  parameters  g  =  (a,  6, 6*,  mi,  7712)^  in  2D 
or  g  =  (ai,  a2, 03,  Vi,  V2,  V2,  s,  mi,  m2,  in  3D;  Another  is  the  local  shape  param- 
eters, that  is,  snake  position  p  =  {pf  }-^i  with  M  being  the  number  of  sampling 
points  on  the  snake.  We  collect  the  global  and  local  shape  parameters  in  a  vector 
q  =  (g^,P^)^-  In  the  model  fitting,  we  need  to  find  the  optimum  configuration  of 
q  so  that  the  total  energy  in  Eq.  (4.9)  is  minimized.  The  estimation  of  both  global 
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and  local  parameters  in  the  model  fitting  makes  the  minimization  problem  highly 
nonlinear. 

5.1    Related  Work 

A  primary  approach  to  solve  energy  minimization  problems  is  the  calculus  of 
variations.  In  this  approach,  finite  diff'erence  or  finite  element  method  is  used  to 
solve  the  Euler-Lagrange  equations  derived  from  the  energy  functional.  This  method 
can  be  regarded  as  gradient  descent  along  the  potential  energy  surface  of  the  model. 
Vemuri  et  al.  [35]  used  a  gradient  descent  (GD)  method  that  iteratively  updates  the 
shape  parameter  q, 

^''^  =   -  (5.1) 

where  q''^^  and  q*^  are  the  status  of  the  shape  parameter  vector  q  at  iteration  k  and 
k  +  1  respectively.  A  is  a  diagonal  matrix  containing  step  sizes  and  consequently 
effecting  the  speed  of  convergence.  Partial  derivatives  of  the  total  energy  ^^^^  can 
be  represented  by  the  sum  of  the  partial  derivatives  of  the  prior  model  and  data 
energies  defined  in  Chapter  4.  While  the  iterative  equation  (5.1)  may  eventually 
converge  to  the  correct  estimate,  in  practice  it  may  be  unacceptably  slow.  In  addition, 
gradient  descent  does  not  guarantee  finding  the  global  minimum  if  the  energy  surface 
is  not  convex  and  the  initial  contour  model  is  far  from  the  target. 

Amini  et  al.  [2]  use  dynamic  programming  (DP)  to  carry  out  a  more  extensive 
search  for  a  global  minima.  Although  DP  provides  numerical  stability,  the  stability 
is  achieved  at  a  high  cost  in  computational  complexity.   Williams  and  Shah  [85] 
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proposed  an  alternative  greedy  algorithm  to  DP  which  drastically  cuts  numerical 
costs,  however,  it  does  not  guarantee  numerical  stability. 

Blake  and  Zisserman  [6]  propose  an  graduated  non-convexity  algorithm,  which 
is  based  on  iterative  approximation  of  the  energy  functional  by  a  convex  function.  It 
can  bridge  low  ridges  in  the  potential  field  and  eliminate  termination  in  local  valleys. 

Poon  et  al.  [57]  and  Grzeszczuk  and  Levin  [25]  minimize  the  energy  of  ac- 
tive models  using  simulated  annealing  (SA)  method  which  is  known  to  give  global 
solutions  and  allows  the  incorporation  of  non-differentiable  constraints,  however,  de- 
terministic algorithms  usually  preferred  due  to  their  faster  convergence. 

5.2    Our  Fast  Numerical  Algorithms 

In  the  following  sections,  we  present  two  algorithms  for  solving  the  energy  min- 
imizing problem  efficiently.  The  first  one  involves  a  Preconditioned  Nonlinear  Con- 
jugate Gradient  (PNCG)  method  for  estimating  the  global  and  local  parameters  of 
the  model.  The  preconditioner  we  use  for  this  case  is  the  diagonal  Hessian  for  the 
global  as  well  as  local  parameters.  In  the  second  numerical  scheme,  we  first  solve 
for  the  global  parameters  using  the  Levenberg-Marquardt  (LM)  method  in  the  outer 
loop  and  employ  the  alternating  direction  implicit  (ADI)  scheme  in  the  inner  loop 
to  solve  for  the  local  parameters.  The  first  scheme  is  less  sophisticated  and  easier  to 
program,  the  second  scheme  is  more  sophisticated  and  relatively  difficult  to  program. 
The  next  sections  will  be  devoted  to  the  presentation  of  the  details  of  each  of  these 
schemes. 
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5.2.1    Preconditioned  Nonlinear  Conjugate  Gradient  Algorithm 

In  general,  the  Conjugate  Gradient  (CG)  algorithm  is  well  suited  for  solving 
quadratic  optimization  problems  whereas,  the  nonlinear  version  can  be  used  for  find- 
ing the  locally  optimal  solution  of  a  nonlinear  non-convex  function  [66].  In  this  work, 
we  develop  a  Preconditioned  Nonlinear  Conjugate  Gradient  (PNCG)  algorithm  and 
apply  it  to  the  energy  minimization  discussed  in  the  previous  section. 

Given  any  n-variate  continuous  function  /(x)  =  f{xi,X2,...,Xn)  with  gradient 
/'(x)  =  [af^/(x)  5f;/(x)  ...  af;;/(x)]^,  we  can  use  Nonlinear  Conjugate  Gradient 
(NCG)  algorithm  to  minimize  it.  The  residual  r,  =  — /'(x,)  indicates  how  far  we  are 
from  the  correct  value  that  minimizes  /(x),  where  i  denotes  the  number  of  iterations. 

We  use  do,di,  ,dn_i  to  denote  a  set  of  orthogonal  search  directions,  ai  and  ^i 

the  step  sizes,  and  we  choose  a  preconditioner  M  that  approximates  /"(x)  and  has 
the  property  that  M~^r  is  easy  to  compute  for  any  vector  r  [36].  Let  imax  be  the 
maximum  number  or  iterations  and  e  be  the  residue  threshold  for  the  function  /(x), 
the  outline  of  the  Preconditioned  Nonlinear  Conjugate  Gradient  method  is  as  follows: 

1.  Let2  =  0,ro  =  -/'(xo). 

2.  Calculate  a  preconditioner  Mq  ~  /"(xq). 

3.  do  =  Mo'ro. 

4.  Find  ai  that  minimize  /(xj  -I-  ttidj). 

5.  Xi+i  =  Xi  +  Oidi. 

6.  Calculate  the  preconditioner  Mi+i  ^  /"(xi+i). 
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7.  ri+i  =  -/'(xi+i). 


=  max  (  ^     rTM.-^ri   ' 


9.  di+i    =   Ti+i   +  A+ldi. 

10.  While  [i  <  imax  or  \f  \  <  e)  i  =  i+1,  go  to  step  4.  :■] 

In  step  8,  using  the  first  expression  for  leads  to  the  Fletcher- Reeves  method, 
and  using  the  second  formula  leads  to  the  Polak-Ribiere  method.  In  step  4,  the 
function  /(x  +  ad)  is  approximately  minimized  by  setting  a  =  —^rjw^-  Where, 
/"(x)  is  the  Hessian  matrix  oi  the  continuous  function  /(x). 


H 


OXiOXi 

£7X2  ax  1 


dxjidxi  dxnO: 
In  our  case,  it  can  be  written  as 


OX1OX2 
OX2OX2 


'X2 


dxi 


d. 


^X 


X2OXJ1 


(5.2) 


H 


d'^E 

d'^E 

dqidqi 

dqidq2 

dqidqn 

d'E 

d^E 

d^E 

dq2dqi 

dq2dq2 

dq2dqn 

d'E 

d'^E 

d^E 

(5.3) 
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Since  it  is  expensive  to  compute  the  Hessian  at  each  iteration,  to  perform  an  exact  line 
search  without  computing  /"(x),  the  Secant  method  is  used  leading  to  the  following 
approximation 

 [f  (x)]^d 

^[/'(x  +  ad)Fd-[/'(x)rd  ^^-^^ 

Typically,  we  will  choose  an  arbitrary  small  nonzero  number  for  a  on  the  first  Secant 
method  iteration;  on  subsequent  iterations  we  will  choose  x  +  ad  to  be  the  value  of 
x  from  the  previous  Secant  method  iteration.  In  other  words,  if  we  let  a;  denote  the 
value  of  a  calculated  during  Secant  iteration  i,  then  (Ti+i  =  — aj. 

In  NCG,  r(z)  =  -/'(x,)  and  Gram-Schmidt  conjugation  of  the  residuals  is  used 
to  compute  the  search  directions.  As  for  computing  /?,  two  choices  are  given  in  the 
above  algorithm.  The  Fletcher-Reeves  technique  converges  for  a  good  initial  guess 
close  to  the  desired  optimum.  Polak-Ribiere  scheme  often  converges  much  faster 
however,  can  cycle  infinitely  without  converging  in  rare  cases.  Convergence  in  this 
case  is  guaranteed  by  choosing  P  =  max{j3^^,  0).  In  our  implementation,  we  use  the 
Polak-Ribere  scheme. 

When  it's  too  expensive  to  compute  the  full  Hessian  matrix,  it  is  reasonable  to 
use  the  diagonal  of  the  Hessian  as  a  preconditioner.  However,  if  the  initial  guess  is 
too  far  from  a  local  optimum,  the  diagonal  elements  may  not  all  be  positive.  Since  a 
preconditioner  should  be  positive  definite,  non-positive  diagonal  entries  will  violate 
this  property  and  hence  can  not  be  allowed.  A  conservative  approach  in  this  case  is 
to  use  the  identity  matrix  as  a  preconditioner  (M  =  I),  that  is,  not  to  precondition. 
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In  the  following  section,  we  describe  a  semi-numerical  approximation  for  computing 
the  Hessian  of  the  energy  function. 

5.2.2    Computation  of  the  Approximate  Hessian 

In  this  section,  we  will  introduce  a  semi- numerical  method  for  approximating 
the  Hessian  of  the  energy  function  by  dissecting  the  structure  of  the  Hessian  matrix. 
In  the  Hessian  matrix  formula,  it  can  be  seen  that  the  elements  in  one  row  of  Hessian 
matrix,  say  row  i  can  be  obtained  by  taking  the  partial  derivative  with  respect  to  gi 
of  the  vector  (||)^,  where  gi  is  the  ith  entry  of  the  vector  g.  Therefore,  the  ith  row 
of  the  Hessian  has  the  form: 


where  n  is  total  number  of  global  parameters  representing  the  snake  pedal  model.  In 
3D,  the  sub-block  corresponding  to  the  snake  pedal  global  aspect  ratio  parameters 
gs  =  [ai  02  aa]^  has  the  following  form: 


d     dE   dE  dE 


d_  r 

%  [  5g 


(5.5) 


%  L        dg2  "'  dgn 


d     dE   dE  dE 


da\  I  da\   da2  da^ 


dxi  dx2  dxs 


dE  dE  dE 
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d'^E  d'^E  d^E 
dx\  dxl 


dx] 
dai 


dx2 
dal 


dx^ 
da: 


ai  J 


dx-[ 


dxi 
~5a2 


8x2  3x2 


dx] 

0X2 


8x2,  ^3:3 

da[    da2  dai 


,  (5.6) 


where  xi,X2,X3  are  the  world  coordinates  in  which  the  model  is  described. 
Similarly,  the  part  corresponding  to  ge  =  [vl  v2  v3  s]^  is, 


d 

dvi 


dE  dE  dE  dE 
dvi   Bv^  5u3  ~Us 


dE  dE  8E 
dxi  dx2  dxs 
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d'^E  d'^E  d'^E 
dxf  dx'2  dxi 


dx-[ 
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dx2 
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3x2, 
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dx\  dx\  dx-[  dx-[ 
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3vi  3v2  3vl  3s 
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We  can  get 


dE  dE  dE  3E 
3vi  3v2  3v3  ~5s 


3 


3E  3E  3E  3E 
3vi  3v2  3v3  ~ds^ 


,  by  carrying 


out  the  similar  procedure.  Since  we  require  only  the  diagonal  entries  of  the  Hessian, 
only  the  following  elements  need  to  be  computed,  ^(|^),  ^(|§)>  ^(f^) 
^(^).  Observe  that  the  computational  complexity  will  be  tremendously  decreased 
if  only  those  diagonal  elements  are  required  to  be  computed.   From  Eq.    (4.6)  in 


Section  4,  we  can  obtain  the  analytical  form  for  computing 


3E  3E  3E 
3xi  3x2  3x3 


and 
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However,  there  is  no  analytical  form  for  ^  and  Qg.Qg.-: 
use  the  forward  difference  to  compute  these  partial  derivatives  numerically.  Our 
experiments  show  that  this  semi-numerical  method  is  very  efficient  and  convenient 
for  the  computation  of  the  diagonal  of  Hessian  matrix. 

5.2.3    LM  Method  for  Global  Parameter  Estimation 

Although  the  Preconditioned  Nonlinear  Conjugate  Gradient  (PNCG)  method 
performs  quite  well  in  the  global  parameter  minimization,  however,  we  found  that 
the  Levenberg-Marquardt  (LM)  method  has  larger  convergence  range  and  converges 
faster  than  the  pure  PNCG  algorithm.  Therefore,  in  our  second  model  fitting  scheme, 
we  utilize  the  LM  algorithm  to  estimate  the  global  shape  parameters.  In  the  local 
shape  parameter  estimation,  by  representing  the  stiffness  matrix  in  a  Schur  com- 
plement formula  [53],  we  use  the  alternating  direction  implicit  method  in  the  third 
model  fitting  scheme  to  achieve  a  0{N)  convergence  rate  in  the  local  parameter  es- 
timation, where  N  is  the  number  of  nodes  in  the  finite  element  discretization  of  the 
snake  positions.  Fig.  (5.1)  shows  the  framework  of  the  second  and  third  model  fitting 
schemes  [37,  85].  We  now  first  introduce  the  Levenberg-Marquardt  method  in  the 
global  shape  parameter  minimization. 

The  LM  method  essentially  combines  the  Newton  method  and  the  gradient 
descent  (GD)  method  in  solving  the  nonlinear  function  minimization  problem.  When 
the  current  solution  is  far  from  the  true  solution,  the  LM  method  adopts  the  gradient 
descent  (GD)  algorithm.  On  the  other  hand,  when  the  current  solution  is  close  to 


d'^E  d^E  d^E 
'dxf  'd^  'dx^ 
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While  (  Energy  >  e  ) 


One  iteration  of 


Levenberg  -  Mai'quardt  Mcthcxl 
(Ref :  Numerical  Recipes) 


solving  global  parameters 


Solve  local  parameters 


Aiternating  Directional  Implicit  Method 


completely 


(Ref  :  Lai  &  Vemuri  "  97  LAA) 


(a) 


Figure  5.1.  Big  picture  of  numerical  schemes  for  model  fitting. 


the  true  solution,  the  nonlinear  function  can  be  approximated  by  a  quadratic  form 
and  the  LM  method  switches  to  the  Newton  method. 

We  now  adopt  the  LM  method  to  our  model  fitting  problem.  Assume  that 
we  only  need  to  fit  the  snake  pedal  model  to  a  set  of  3D  data  points,  then  the 
external  energy  used  here  becomes  the  spring-based  deformation  energy  Ed(D,  g,  p) 
introduced  in  Section  4.  More  specifically,  to  fit  the  snake  pedal  to  a  given  data  set 
of  m  points  {D;}^^  in  3D,  we  minimize  the  sum  of  the  squared  distances  from  each 
data  point  Dj  to  the  corresponding  closest  point  Xi  on  the  pedal  curve/surface.  Let 
h  be  the  function  that  returns  the  closest  point  on  the  snake  pedal  from  a  given  data 
point  D,,  that  is,  x^  =  hi(g,pi,ij).  Let  {pi}^^  be  a  set  of  M  discrete  points  along 
the  snake  p{t),  {xi}fii  be  the  corresponding  points  on  the  snake  pedal  curve.  Let 
D  =  {Df       and  p  =  {pj}fii,  the  quantity  to  be  minimized  for  the  model  fitting 
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can  be  written  as 


m 

ED(D,g,p)  =  -;^  c(D,-x,)^ 


(5.7) 


where  Xj  =  hi(g,  p,,  tj),  and  the  magnitude  of  parameter  c  is  related  to  the  uncer- 
tainty of  data  measurements. 

In  LM,  the  gradient  of  ED(D,g,p)  with  respect  to  generator  parameters  g  = 
{91,92, --V  =  {aua2,a3,vi,V2,V3,s, mi, 7712,1713)'^  in  3D,  and  g  =  {91,92, ---V  = 
(ai,  02,  9,  nil,  1712)^  in  2D  has  components: 


dE 


D 


dgk 


c(Di-hi(g,pi,i,)) 


1=1 


dhi{g,pi,ti) 
dgk 


(5.8) 


If  we  assume  that  there  is  a  one-to-one  correspondence  between  the  points  on  the 
generator  and  the  pedal  snake  then  the  partial  derivative         can  be  written  as: 


  -  c  (Dk  -  hk(g,Pk,4))- 


dgk 


d9k 


(5.9) 


Taking  an  additional  partial  derivative  gives  Qg^^  which  is  needed  in  the  LM 
iteration.  For  the  derivative  of  the  energy  required  in  the  LM  algorithm,  we  use  the 
chain  rule  to  get  =  (%f  )(§|)-  More  specifically,  if  we  denote  each  component 
hj  of  h  as  hi  =  {hii,hi2Y ,  then  in  2D 


as 


Dii  -  hii    Di2  -  h 


i2 


dhj]  dhu  dhji  dh^  dhu 

da  db  dU  drrli  81712 

dhn  dhg  dhg  dhi2  dh,,2 

da  db  36  dvrii  dml  . 


(5.10) 
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If  we  define 


Pk   =   -^  >  (^kl=  ^       ,  5.11 


using  the  Newton  method  to  minimize  the  energy  function  gives 


J2^ki5gk  =  Pk-  (5.12) 


On  the  other  hand,  using  the  gradient  descent  method  gives 


Sgi  =  const,  x  /3i  =  (5.13) 


where  we  choose  the  constant  to  be  ^  in  Eq.  (5.13).  This  choice  of  step  size  in  Eq. 
(5.13)  leads  to  a  preconditioned  gradient  descent  method  with  a  diagonal  Hessian  as 
the  preconditioner.  Furthermore,  we  define 

4  =  «J-^(1  +  ^)'  (5-14) 

"jfc  =  "jfc  (Jt^^),  (5.15) 
and  combining  Eq.  (5.12)  and  (5.13)  yields 


^oi'f.i6gi  =  pk,  (5.16) 


When  A  is  very  large,  Eq.  (5.16)  becomes  identical  to  (5.13);  and  as  A  approaches 
to  zero,  Eq.  (5.16)  becomes  the  same  as  Eq.  (5.12). 
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From  above  discussion,  given  an  initial  guess  for  the  set  of  parameters  g,  the 
LM  method  for  the  global  parameter  estimation  is  described  as  follows: 

1.  Compute  Ed(D, g, p), 

2.  Pick  a  suitable  value  for  A,  for  example,  A  =  0.001, 

3.  (t)  Solve  the  linear  Eq.  (5.16)  for  Sg  and  evaluate  ED(D,g  +  5g,p), 

4.  If  Ed(D,  g  +  6g,  p)  >  Ed(D,  g,  p),  increase  A  by  a  factor  of  10  (or  any  other 
substantial  factor)  and  go  back  to  (f), 

5.  If  ED(D,g  +  (^g,p)  <  ED(D,p,q),  decrease  A  by  a  factor  of  10,  update  the 
trial  solution  g  < —  g  +  Sg,  and  go  back  to  (f),  until  the  energy  Ed(D,  g,  p)  is 
less  than  a  predefined  threshold  . 

After  the  estimation  of  the  global  parameters,  we  solve  the  local  shape  parame- 
ters in  the  inner  loop  of  the  parameter  estimation  using  the  Preconditioned  Conjugate 
Gradient  method  or  the  alternating  directional  implicit  method. 

5.2.4    API  Method  for  Local  Parameter  Estimation 

When  solving  for  the  local  parameters  in  the  inner  loop  for  fixed  global  pa- 
rameters, the  quantity  to  be  minimized  has  a  quadratic  form  and  its  minimization 
leads  to  solving  a  linear  system,  namely,  Kx  =  b,  where  K  is  the  stiffness  matrix, 
X  is  the  local  parameter  vector,  and  b  =  is  the  external  force.  Note  that  we 

use  X  to  represent  the  local  shape  parameter  vector  and  b  to  represent  the  external 
force,  (they  were  represented  by  p  and  f  in  previous  chapters.)  which  is  the  standard 
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method  to  represent  unknowns  and  right-hand  side  in  linear  algebra.  Because  of  the 
special  structure  of  the  matrix  K  in  this  linear  system,  we  can  use  the  alternating 
direction  implicit  (ADI)  method  in  the  local  shape  parameter  estimation.  We  can 
prove  that  the  ADI  takes  only  0{N)  time  to  converge  with  soft  data  constraints 
where  N  is  the  number  of  nodes  in  the  rectangular  finite  element  discretization  of 
the  snake  in  our  model.  Note  that  in  soft  data  constraints,  the  model  is  not  required 
to  interpolate  the  data  but  is  expected  to  approximate  it. 

We  now  present  the  structure  of  the  stiffness  matrix  K.  In  the  snake  pedal 
model,  the  controlling  snake  p  is  a  controlled  continuity  spline  under  the  influence  of 
image  forces  or  other  external  constraint  forces.  In  our  model,  the  spline  deformation 
energy  imposed  on  the  snake  is  a  membrane  plus  thin  plate  energy  given  in  continuous 
form  by 


where  Wi  and  W2  are  the  tension  and  rigidity  factors  respectively. 

The  external  energy  used  can  be  the  same  as  the  energy  Ed(D,  q,  p)  discussed 
in  Section  5.2.3.  In  this  case,  the  external  force  impose  on  the  snake  is  given  by 


^,(P)  =  /».((|)^  +  (|)^)  +  Mi^?  ^  2i^/  +  (0)'^)  <lu  (5.17) 


b  = 


c  [Da  -  ha  A2  -  ^1,2]  ■ 


dhi 


(5.18) 
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South  Pole 


Total  number  of  nodes 

n  X  n  +  2 
Denote  :  N  =  n  X  n 


North  Pole 


(a) 


Figure  5.2.  Model  tessellation  in  the  model  coordinates. 


where 


dhi2 

2   I  2 


"1  ■ 

i  —       r>    17 

n(  +  n2 


(5.19) 


with  ni,n2  defined  as  before. 

Notice  that  dEo/dhi  is  the  external  force  in  a  standard  snake,  but  in  our  case, 
the  external  force  is  dEo/dpi,  which  is  related  to  the  force  in  a  standard  snake  via 
a  Jacobian  matrix  dhi/dpi  in  Eq.  (5.19). 

We  then  discretize  the  model  in  model  coordinates  using  finite  elements.  Ge- 
ometrically, the  snake  pedal  surfaces  are  closed  surfaces  in  space  whose  intrinsic 
coordinates  u  =  {u,  v)  are  defined  on  a  domain  1].  The  positions  of  points  on 
the  model  in  an  inertial  frame  of  reference  are  given  by  a  vector-valued  function 
pedal(u)  =  (a;i(u),Z2(u),X3(u)).  The  tessellation  of  u  =  {u,v)  in  intrinsic  coordi- 
nates is  shown  in  Fig.  5.2.  We  use  a  combination  of  quadrilateral  and  triangular 
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elements  in  the  discretization.  In  the  inner  quadrilateral  mesh,  we  have  {nx  n)  =  N 
nodes  and  2  additional  nodes  are  used  in  representing  the  north  and  south  pole  lead- 
ing to  a  total  of  {N  +  2)  nodes.  Due  to  the  presence  of  the  north  and  south  poles, 
after  the  discretization  of  the  model  in  intrinsic  coordinates  u  using  finite  elements, 
the  stiffness  matrix  K  E  i?(^+2)(Ar+2)  -j^  ^^j.  problem  can  be  written  as  a  (2  x  2)  block 
matrix: 


K  = 


Kii  Ki2 
K21  K22 


(5.20) 


where  Kn  G  fi^^^  contains  coefficients  which  corresponds  to  the  inner  quadrilateral 
mesh  and  K22  £  R^^^  is  a  diagonal  block  matrix  containing  the  coefficients  corre- 
sponding to  the  north  and  south  poles,  and  K12  =  Kji  contains  the  coefficients  of 
interaction  between  the  poles  and  other  nodes. 

To  solve  the  linear  system  Kx  =  b  for  the  local  parameter  estimation,  we  employ 
the  following  Schur  complement  formula  [53] 


K-^b 


Kr/  +  Kr/Ki2s-iK2iKr/  -Kr/Ki2s-i 


(5.21) 


where 


S  —  K22  —  K21K1/K12. 


(5.22) 
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From  Eq.  (5.21),  it  is  evident  that  to  obtain  the  solution  of  K~^b,  we  need  the 
solution  of  several  K^i  bs.  Further  analysis  (see  Appendix  A)  shows  that  the  com- 
ponent Kii  of  the  stiffness  matrix  K  can  be  expressed  as  the  sum  of  two  Kronecker 
products: 

Kii   =   u;i(A(8)B  +  B(8)  A) +  i(;2A(8)  A 

=  A(8)C  +  C®A,  (5.23) 

where  C  =  +  {w2/2)A  with  A  and  B  being  (n  x  n)  tridiagonal  and  symmetric 
positive  definite  (SPD)  matrices  of  the  form 


1  -1 


2  1 


1      2  -1 


1     4  1 


A  = 


-1      2  -1 


1     4  1 


-1 


1  2 


(5.24) 

The  Kronecker  tensor  product  or  Kronecker  product  of  C  and  A,  denoted 
by  C  (8)  A,  is  the  matrix  of  order  N  x  N  (where  N  =  n  x  n)  and  given  by 
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ciiA  C12A 


ClnA 


C21A  C22A 


C2nA 


C0A  = 


(5.25) 


CjiiA.  CJ12A.  ...  CfijjA 
When  computing  with  Kronecker  products,  it  is  computationally  convenient  for 
us  to  represent  vectors  using  matrices.  When  solving  (C  (E>  A)x,  we  represent  the 
vector  X  of  order  •n?  by  the  matrix  X  =  {xki}  of  order  n  x  n  defined  by 


It  can  be  shown  that  this  representation  can  lead  to  efficient  procedures  for 
computing  (C  (g)  A)x  and  solving  (C  (g)  A)x  =  b,  respectively. 

LEMMA  1 

Let  C  =  {cki},  A  =  {cfc/}  and  X  =  {xki}  be  matrices  of  order  n  x  n,  then  the 
n  X  n  matrix  representation  of     x  1  vector  (C  (g)  A)x  is  given  by 


Xkl  —  Xk+n{l-l)- 


(5.26) 


(CO  A)x 


(C(AX)T)T. 


Using  this  representation,  the  linear  system  Knx  =  b  can  be  written  as: 


CXA^  +  AXC^  =  D, 


(5.27) 


.  ■  68 

where  D,X  are  the  n  x  n  representation  of  the  n'^  vector  b  and  x.  Since  A  and  C 
are  also  SPD  matrices,  Eq.  (5.27)  can  simply  be  written  as 

CXA  +  AXC  =  D.  (5.28) 

After  some  simple  linear  algebra  operations,  it  follows  that  above  equation  is  equiv- 
alent to 

A'X  +  XA'  =  D',  (5.29) 

with  A'  =  A-^C  and  D'  =  A'^DA-^ 

Eq.  (5.29)  is  the  standard  Lyapunov  matrix  equation  [43].  The  ADI  method  can 
be  applied  to  solve  this  matrix  equation  and  has  the  following  steps  in  each  iteration 
J  =  1,2,...,  J, 

(A' +  p,I)X^._i   =  D'-X,_i(A'-p,I)  (5.30) 
X,(A'+p,I)   =  D'-(A'-p,I)X^_,.  (5.31) 

For  each  iteration  in  ADI,  we  need  to  solve  two  equations,  namely,  Eq.  (5.30)  and 
(5.31),  and  therefore  X:!  is  introduced  as  an  intermediate  variable.  Since  A'  is 
a  SPD  matrix,  we  can  use  Cholesky  decomposition  for  the  tridiagonal  SPD  matrix 
(A'+Pjl)  to  compute  the  update  solution.  This  only  takes  0{N)  time  per  iteration. 
The  parameters  pj  are  called  ADI  parameters  and  J  is  the  number  of  iterations 
needed.  They  are  chosen  according  to  the  method  described  by  Lai  and  Vemuri  [36]. 
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In  order  to  fit  the  snake  pedal  curve  to  salient  features  in  an  image,  we  use  the 
following  equation  as  the  energy  to  be  minimized: 


E  =  -||(G,*V/)|r  =  -||(G,*V/(x,2/)|p.  (5.32) 

■        ■  , 

Where  (x,  y)  is  on  the  snake  pedal  curve,  is  a  Gaussian  with  standard  deriva- 
tion (7,  and  *  denotes  the  convolution  operator.  More  sophisticated  boundary  finding 
energies  can  be  synthesized  and  we  refer  the  reader  to  Kass  et  al.  [31].  By  mini- 
mizing the  above  energy,  the  model  is  made  to  fit  to  points  of  high  gradient  in  the 
image.  The  procedure  to  minimize  above  energy  is  similar  to  that  of  fitting  of  snake 
pedal  curves/surfaces  to  data  points. 


CHAPTER  6 
HYBRID  GEOMETRIC  ACTIVE  MODELS 


Since  the  inception  of  snakes  in  the  vision/graphics  community  by  Kass  et  al. 
[31]  and  Terzopoulos  et  al.  [76],  these  elastically  deformable  contours/surfaces  have 
been  widely  used  for  a  variety  of  applications,  such  as  edge  detection,  segmentation, 
motion  tracking,  and  shape  modeling.  The  classical  approach  to  object  shape  re- 
covery using  the  snakes  is  based  on  deforming  an  initial  configuration  of  the  snake 
Vo  towards  the  boundary  of  the  object  to  be  detected.  As  described  in  Chapter  2, 
snakes  are  energy-minimizing  splines,  their  deformation  under  constraints  is  achieved 
by  minimizing  a  functional,  that  can  be  regarded  as  the  bending  energy  stored  in  a 
thin  flexible  beam/rod  or  a  stretchable  string  subject  to  loading. 

There  are  several  problems  associated  with  this  approach,  such  as  initializations, 
convergence  to  the  local  minimum,  and  the  specification  of  the  physical  or  elasticity 
parameters.  Moreover,  this  energy  model  requires  that  the  topology  of  the  shape  to  be 
estimated  be  known  a  prior.  However,  in  many  applications,  such  as  motion  tracking, 
it  is  not  always  possible  to  specify  the  topology  of  an  object  prior  to  its  recovery. 
For  example,  during  the  cell  mitosis,  the  close  contours  of  cell  boundary  change 
connectivity  or  split,  therefore  undergoing  a  topological  transformation.  Hence  it  is 
necessary  to  develop  a  model  which  has  the  ability  to  split  and  merge  freely  during  the 
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shape  recovery  without  making  the  a  prior  assumption  about  the  object's  topology 
in  these  applications. 

Recently,  Caselles  et  al.  [9]  and  Malladi  et  al.  [45]  proposed  a  novel  geometric 
model  of  active  contours  dubbed  geodesic/ geometric  active  contours  or  geodesic/ geometric 
snakes.  These  models  are  based  on  the  the  theory  of  curve  evolution  and  geometric 
flows.  Automatic  changes  in  topology  can  be  handled  in  a  natural  way  in  this  model- 
ing technique  when  we  implement  the  curve  evolution  using  the  level-set  embedding 
schemes.  In  essence,  this  model  is  and  implicit  surface-evolution  model  where  the 
surface  propagates  with  a  velocity  which  is  directly  related  to  the  curvature  of  the 
evolving  curve  which  is  embedded  in  the  surface.  Unlike  the  classical  energy  models, 
these  models  are  not  the  result  of  minimizing  an  energy  functional.  In  this  frame- 
work, boundary  detection  can  be  considered  equivalent  to  finding  a  curve  of  minimal 
weighted  length.  More  recently,  several  researchers  independently  demonstrated  the 
unity  of  the  curve  evolution  approaches  and  the  classical  energy  minimization  meth- 
ods [10,  32,  64]. 

The  challenge  in  the  geometric  active  contour  evolution  is  to  devise  numerical 
schemes  for  the  equations  of  the  evolving  curve,  which  will  accurately  approximate 
solution  that  can  be  very  unstable  due  to  formation  of  shocks  or  discontinuities.  Os- 
her  and  Sethian  [52]  developed  special  numerical  techniques  based  on  upwind  finite 
difference  schemes  leading  to  a  very  stable  solution.  They  introduced  level-set  tech- 
niques to  seamlessly  handle  any  topological  changes  of  the  curve/surface  that  may 
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occur  during  evolution.  The  equation  of  motion  for  the  embedded  curve  is  remi- 
niscent of  an  initial  value  "Hamilton-Jacobi"  equation  with  a  parabolic  right-hand 
side  and  is  closely  related  to  a  viscous  hyperbolic  conservation  law.  For  details  on 
the  theory  of  curve/surface  evolution  and  its  level-set  implementation,  we  refer  the 
reader  to  the  relative  papers  [9,  10,  33,  46,  51,  45,  65,  74]. 

6.1    Overview  of  Hybrid  Geometric  Active  Models 

In  many  computer  vision  applications,  characterizing  the  global  shape  of  an  ob- 
ject is  very  useful,  for  example,  in  object  recognition.  Traditional  geometric  active 
contour/surface  models  do  not  possess  the  capability  to  characterize  the  global  shape 
of  an  object.  They  are  however  useful  in  characterizing  the  local  shape  details  in  an 
object.  In  this  section,  we  introduce  a  novel  concept  of  a  global/core  model  into  the 
PDE-based  curve  evolution  framework  by  embedding  the  snake  pedal  model  into  a 
level-set  framework.  Instead  of  characterizing  an  object  boundary  by  the  position  of 
every  point  on  the  boundary,  our  proposed  model,  referred  to  as  the  hybrid  geometric 
active  model,  now  describes  an  object  as  a  combination  of  a  global/core  shape,  for 
example,  circle,  super-ellipse,  and  so  on.  and  a  variable  offset  defined  with  respect 
to  the  global  shape.  The  variable  offsets  are  controlled  by  this  global  shape  and  an 
evolving  curve — the  controlling  snake  in  the  snake  pedal  model.  Augmentation  of 
the  global  shape/core  in  the  model  representation  is  very  useful  in  object  recognition 
and  image  indexing  applications.  In  these  applications,  the  global  description,  such 
as  center,  size,  and  orientation,  of  a  group  of  objects  or  structures  needs  to  be  re- 
trieved initially  which  is  then  followed  by  a  retrieval  of  the  local  description  of  each 
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object  (if  necessary).  For  the  model  to  recover  the  object  boundary,  we  introduce 
a  robust  and  efficient  numerical  method  which  consist  of  a  global  plus  local  shape 
estimation  technique  in  the  model  fitting.  We  use  the  Levenberg-Marquardt  (LM) 
described  in  Chapter  5  for  estimating  the  global  shape  and  a  combination  of  up- 
wind and  minmod  finite  difference  schemes  in  a  level-set  framework  for  estimating 
the  local  shape.  The  LM  method  has  a  broad  convergence  range,  converges  to  the 
local  optimum  very  quickly,  and  therefore  is  well  suited  for  the  global  shape  recov- 
ery. The  hybrid  geometric  active  contour/surface  model  retains  all  the  advantages  of 
traditional  geometric  active  models  (for  example,  topology  change,  ability  to  model 
complex  geometries  and  amenability  to  stable  numerical  implementation)  and  has 
the  added  ability/ advantage  of  being  able  to  compactly  represent  global  shape  of  an 
object. 

We  first  present  the  theoretical  foundation  of  the  traditional  curve  evolution  in 
Section  6.2,  followed  by  the  level-set  implementation  framework  of  the  PDE-based 
curve  evolution  in  Section  6.3.  We  then  introduce  our  hybrid  geometric  active  models 
in  Section  6.4.  The  numerical  issues  of  the  model  fitting  process  will  be  discussed  in 
Chapter  7,  followed  by  the  implementation  results  in  Chapter  8. 

6.2    Curve  Evolution  Theory 

The  mathematical  foundation  of  the  geodesic  or  geometric  active  contour  model 
is  based  on  the  Euclidean  curve  shortening  theory,  which  guides  the  curve  evolution 
in  Euclidean  space  [15].  Let  p  =  {px.PyY  denote  an  arbitrary  point  on  the  geomet- 
ric active  contour,  the  family  of  planar  curves  flow  according  to  the  geometric  heat 
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equation  or  curve  shortening  equation 


where  k  is  the  curvature  and  N  is  the  inward  unit  normal.  When  the  curve  evolves 
according  to  (6.1),  the  Euclidean  perimeter  shrinks  as  quickly  as  possible.  We  demon- 
strate this  concept  in  the  following. 

Let  7^(0)  be  a  smooth,  closed  initial  curve  in  Euclidean  plane  3fJ^,  let  V{t)  be 
the  one-parameter  family  of  curves  generated  by  moving  'P(O).  Let  p{s\t)  be  the 
position  vector  which  parameterizes  V{t)  by  s',  where  0  <  s'  <  1  and  t  is  the  time 
factor.  For  closed  curves,  we  assume  that  p(0,  t)  =  p{l,t)  for  any  time  t  and  a  similar 
condition  is  posed  on  the  first  derivatives.  We  define  the  length  functional: 


(6.2) 


where  {|.||  denotes  the  standard  2-norm.  Taking  the  first  variation  [19]  and  using 


integration  by  parts  in  Eq.  (6.2),  we  obtain 


dp  d'^p 


(6.3) 


as'" 


1  d 


-  dp  - 

ds' 


11011  >  ds', 
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where  <,  >  is  the  inner  product  of  two  vectors.  Note  that  the  arc  length  ds  can  be 
expressed  by 


using  the  definition  of  curvature  k  and  normal  N  [7,  8,  49],  the  integral  in  (6.3) 
becomes 


Setting  equation  (6.5)  to  zero,  we  obtain  the  direction  in  which  C{t)  decreases  most 
rapidly  as 


It  has  been  demonstrated  [16,  17,  18,  24,  60]  that  when  simple  closed  curves  evolve 
according  to  (6.1)  without  developing  singularities,  they  will  converge  to  "round" 
points.  High  curvature  regions  will  be  smoothed  out  during  the  evolution.  We  can 
easily  understand  this  phenomenon  by  further  examining  Eq.  (6.1).  In  (6.1),  the 
planar  curve  evolves  in  the  direction  of  its  inward  normal,  and  the  rate  of  evolution 
is  proportional  to  its  local  curvature.  In  locations  where  the  curve  has  a  larger  cur- 
vature, the  rate  of  change  is  large,  thus  the  curve  evolves  faster  in  the  high  curvature 
points.  For  the  same  reason,  the  curve  evolves  relatively  slower  at  the  low  curvature 
points.  Therefore,  the  portions  of  the  curve  in  high  curvature  regions  undergo  more 
changes  toward  the  rest  state,  while  the  portions  in  the  low  curvature  regions  remain 
smooth  until  all  portions  of  curve  reach  the  rest  state — a  circle. 


ds  =  ||— 7II  ds', 


(6.4) 


(6.5) 
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Since  the  curvature  k  determines  the  rate  of  change  in  (6.1),  it  is  regarded  as 
the  speed  function  of  the  curve  in  the  evolution.  More  general  speed  functions  can  be 
allowed  in  the  curve  evolution.  Let  F{k)  denote  the  speed  function,  the  curve  then 
evolves  according  to 

^  =  nm-  (6.6) 

In  general,  F{k)  can  be  split  into  two  terms: 

F{k)  =  FA  +  FG,  (6.7) 

where  the  first  term  Fa  is  referred  to  as  the  advection  term  and  is  independent  of 
the  geometry  of  the  curve.  The  curve  uniformly  expands  or  contracts  with  speed 
Fa  depending  on  its  sign,  the  behavior  of  the  curve  is  analogous  to  that  of  balloon 
model  [12]  with  an  inflation/deflation  force.  The  second  term  Fq  is  referred  to  as 
the  diffusion  term  and  is  dependent  on  the  geometry  of  the  curve,  such  as  the  local 
curvature.  The  diffusion  term  smoothes  out  the  high  curvature  regions  of  the  curve 
and  has  the  same  regularization  effect  on  the  curve  as  the  internal  deformation  energy 
term  in  the  generalized  splines  described  in  Chapter  2. 

When  implementing  the  curve  evolution  equations  such  as  (6.6),  there  are  a 
number  of  numerical  considerations.  One  numerical  approach  is  first  taking  the 
above  Lagrangian  description  of  the  problem,  producing  equations  of  motion  for  the 
position  vector  p{s',t),  then  discretizing  the  parameterization  with  a  set  of  discrete 
markers  lying  on  the  moving  curve.   These  discrete  markers  are  updated  in  time 
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by  approximating  the  spatial  derivatives  in  the  equations  of  motion  and  advancing 
their  positions.  However,  there  are  several  problems  with  this  approach,  as  discussed 
in  Sethian  [61].  First,  small  errors  in  the  computed  marker  positions  are  tremen- 
dously amplified  by  the  curvature  term,  and  calculations  are  unstable  unless  an  ex- 
tremely small  step  is  employed.  Secondly,  singularities  may  develop  (when  Fa  7^  0) 
in  the  propagating  curve  even  with  a  smooth  initial  curve.  A  typical  example  of 
this  phenomenon  is  when  the  smoothing  curvature  (viscous)  term  is  absent,  namely, 
F{k)  =  F4  in  (6.6),  and  an  entropy  condition  must  be  imposed  to  extract  the  correct 
weak  solution.  Thirdly,  it  is  difficult  to  handle  the  topological  changes  if  the  evolv- 
ing curve  breaks  and  merges.  Finally,  significant  bookkeeping  problems  occur  in  the 
extension  of  this  technique  to  three  dimensions. 

6.3    Traditional  Level-Set  Approach 

Osher  and  Sethian  [52,  61]  proposed  an  algorithm  based  on  the  theory  of  vis- 
cosity solutions  to  PDEs  .  The  algorithm  provides  a  stable  numerical  solution  for 
curve/surface  evolution.  In  their  approach,  the  curve  is  first  embedded  in  a  two 
dimensional  surface,  the  equations  of  motion  are  solved  using  finite  diflference  tech- 
niques and  hyperbolic  conservation  laws. 

6.3.1    Embedding  Mechanism 

The  embedding  step  can  be  carried  out  by  representing  the  curve  V{t)  by  the 
zero  level  set  of  a  smooth  and  Lipschitz  continuous  function  ^  :  x  [0,  r)  ->  As 
an  illustration  of  this  approach,  we  consider  the  example  of  an  expanding  circle.  Let 
the  initial  curve  P  at  ii  =  0  be  a  circle  in  the  xy-plane  (Fig.  6.1  (a)).  We  imagine 
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Figure  6.1.  Level  set  formulation  of  equations  of  motion  -  (a)  and  (b)  show  the 
curve  V  and  the  surface  (i){x,y)  at  i  =  0,  (c)  and  (d)  show  the  curve  V  and  the 
corresponding  surface  0(a;,  y)  at  time  t.  (e)  and  (f )  show  the  curve  V  and  surface 
</)(x,  y)  at  time  t  +  Ai.  The  curve  V  has  changed  topology  in  going  from  (c)  to  (e). 
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that  the  circle  is  the  level  set  {0  =  0}  of  an  initial  surface  z  =  (f){x,y,t  =  0)  in  dt^ 
(Fig.  6.1(b)).  We  can  then  match  a  one-parameter  family  of  moving  curves  Vit)  with 
a  one-parameter  family  of  moving  surfaces  in  such  a  way  that  the  level  set  =  0} 
always  yields  the  moving  curve  (Fig.  6.1  (c)  and  (d)). 

6.3.2    Equation  of  Motion 

In  the  general  case,  let  7^(0)  be  a  closed,  non-intersecting,  N  —  1  dimensional 
hypersurface.  Let  0(p,i),  p  6  be  a  scalar  function  such  that  0(p,O)  =  ±d{p), 

where  d{p)  is  the  signed  distance  from  p  to  the  hypersuface  'P(O).  Assume  that  0  is 
negative  in  the  interior  and  positive  in  the  exterior  of  the  zero  level  set. 

Let  7^(0)  be  a  planar  curve,  represented  as  the  zero  level  set  of  (j){p,t),  that  is, 

V{t)  e      :  0(p,t)  =  0,  (6.8) 
By  differentiating  (6.8)  with  respect  to  t  we  obtain: 


V(j>{p,t)-^  +  MP,t)=0.  (6.9) 


Note  that  the  gradient  V(j){p,  t)  is  normal  to  the  N  -  1  dimensional  level  set  passing 
through  p  in  general.  When  A^^  =  3,  for  the  zero  level  set,  the  following  relation  holds: 
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where  N  is  the  curve  normal  as  defined  before.  In  this  equation,  the  left  hand  side  is 
in  terms  of  the  surface  (f),  while  the  right  hand  side  is  in  terms  of  a  property  of  the 
curve  V.  Combining  Eq.  (6.1),  (6.10),  and  (6.9)  yields: 


(j>t  =  F{k)\\V^\\.  (6.11) 

The  curve  V,  evolving  by  Eq.  (6.6),  is  obtained  by  the  taking  the  zero  level  set 
of  the  function  0,  which  evolves  by  (6.11).  This  scheme  is  referred  to  as  Eulerian 
formulation  for  curve  evolution,  because  it  is  written  in  terms  of  a  fixed  coordinate 
system.  Also,  for  an  implicitly  defined  curve,  the  curvature  can  be  computed  as 
Morgan  [49] 


or  explicitly  as 


Eq.  (6.11),  together  with  the  initial  condition,  (f){p,0)  =  ±d{p),  is  the  governing 
equation  for  the  curve  evolution  in  a  level-set  framework.  Since  this  equation  re- 
sembles the  Hamilton-Jacobi  equation,  it  is  referred  to  as  the  "Hamilton-Jacobi" 
formulation  [52]. 

To  attract  the  active  contour  to  the  desired  image  features,  the  speed  function 
in  (6.11)  is  multiplied  by  a  quantity: 


81 


where  /  is  the  grey-scale  image  and  Ga  is  the  Gaussian  smoothing  filter  with  width 
(7,  and  n  =  1  or  2  depends  on  the  application  [44,  83].  The  function  g{x,y)  depends 
on  the  given  image  and  is  used  as  a  "stopping  term."  For  example,  when  the  curve 
is  close  to  an  edge,  g{x,y)  becomes  very  small  and  retards  the  curve  evolution.  A 
modification  of  (6.11)  based  on  the  use  of  Riemanian  metric  was  reported  in  Caselles 
et  al.  and  Kichenassamy  et  al.  [10,  32],  it  involves  the  minimization  of  a  weighted 
Euclidean  length  and  takes  the  form: 

(l>t  =  g{x,  y)  F{k)  II       +  V0  •  V^,  (6.15) 

Note  that  for  g  as  in  (6.14),  Vg  behaves  like  a  "doublet"  near  an  edge.  The  effect 
of  Vg  is  to  attract  the  evolving  contour  as  it  approaches  an  edge,  and  to  push  the 
contour  back  if  it  passes  the  edge.  The  inclusion  of  ■  Vg  term  allows  the  use  of 
larger  inflation  forces,  resulting  in  features  being  extracted  in  relatively  fewer  time 
steps. 

6.3.3    Advantages  of  the  Level-set  Approach 

By  using  the  level  set  approach  discussed  above,  topological  changes  can  be 
handled  naturally.  This  advantage  lies  in  the  fact  that  in  the  curve  evolution,  the 
higher  dimensional  function  0(p,  t)  always  remains  a  function,  while  its  level  surface 
{0  =  0},  which  corresponds  to  the  evolving  curve  p{t),  needs  not  to  be  simply 
connected.  The  moving  curve  can  change  its  topology  or  form  sharp  corners  easily.  As 
shown  in  Fig.  6.1  (f),  (f>{p,t)  remains  a  simple  continuous  saddle-like  function,  while 
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its  level  set  {0  =  0},  the  moving  curve,  splits  into  two  separated  closed  contours,  thus 
undergoes  a  topological  change  from  its  initial  configuration  shown  in  Fig.  6.1  (a). 
Another  advantage  of  the  level  set  approach  lies  in  the  fact  that  the  geometric  and 
differential  properties  of  p{t)  can  be  computed  from  0.  For  example,  the  curvature 
of  p{t)  can  be  computed  in  terms  of  the  derivatives  of  0  expressed  in  Eq.  (6.13). 
In  addition,  the  level  set  approach  can  be  easily  extended  to  higher  dimensions  and 
appropriate  expressions  can  be  obtained  for  the  mean  curvature  and  the  Gaussian 
curvature.  Moreover,  the  governing  equation  (6.11)  is  reminiscent  of  an  initial  value 
Hamilton- Jacobi  equation  with  viscosity,  and  can  be  solved  using  the  stable,  entropy- 
satisfying  finite  difference  schemes  in  accordance  with  the  hyperbolic  conservation 
laws  [40,  52,  62].  , 

6.4    Our  Model:  Hybrid  Geometric  Active  Contours 

As  described  in  the  former  sections,  the  level  set  approach  provides  an  elegant 
way  to  handle  topological  changes  with  active  models,  however,  traditional  PDF- 
based  approaches  do  not  possess  a  mechanism  to  characterize  the  global  shape  of 
an  object.  In  this  section,  we  describe  how  the  level  set  formulation  of  the  curve 
evolution  can  be  applied  to  our  snake  pedal  model  to  realize  topological  changes 
(when  necessary)  as  well  as  capture  the  global  shape  representation  of  an  object. 
6.4.1    Big  Picture 

In  our  approach,  to  build  some  amount  of  regularization  on  the  snake  pedal, 
we  impose  the  regularization  constraint  via  Euclidean  arc-length  minimization  on 
the  snake  pedal.  This  would  lead  to  the  standard  geodesic  active  contour.  Let  us 
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consider  a  snake  pedal  denoted  by  Ve{t),  with  its  position  denoted  by  Pe  =  {pei,Pe2}) 
then  the  standard  curve  evolution  formula  for  the  snake  pedal  can  be  written  as 

^  =  F{k,)N,,  (6.16) 

where  ke  and  Ng  are  the  curvature  and  normal  of  the  snake  pedal  respectively.  We 
would  like  to  know  how  the  "snake"  would  evolve  if  the  snake  pedal  were  evolving  as 
a  function  of  its  local  curvature.  Indeed,  we  will  not  evolve  the  snake  pedal  curves 
directly,  instead,  we  will  first  solve  the  "snake"  position  under  the  constraints  that 
the  arc-length  of  the  snake  pedal  is  minimized,  and  the  pedal  curves  can  then  be 
determined  by  the  pedaling  operation  defined  in  Chapter  3,  given  the  position  of  the 
generator.  The  problem  therefore  can  be  solved  by  the  following  procedure: 

1.  Derive  the  curve  evolution  equation  for  the  "snake"  'P{t)  by  minimizing  the 
arc-length  of  the  snake  pedal  'Pe(t); 

2.  Embed  the  evolving  curve  Vit)  in  a  higher  dimensional  surface  0i  and  formulate 
the  equation  of  motion  for  <pi ; 

3.  Solve  the  equation  of  motion  for  0i  using  proper  numerical  techniques; 

4.  Determine  the  snake  pedal  curve  from  the  evolving  snake  via  the  pedaling 
operation,  given  the  generator. 

We  should  remind  the  reader  that  in  the  above  procedure,  the  snake  pedal 
curve  also  evolves,  but  it  evolves  as  a  standard  geometric  active  contour  embedded 
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P  e  (t)  =  Level  Set  ^^^=0 


;?(t)  =  Level  Set    <t)  =0 


y 


(a) 


Figure  6.2.  Relation  between  V{t),Ve{t),(l)i,  and  (?!)2:  V{t)  and  'Pe(i)  are  level-sets 
of  the  higher  dimensional  surface  and  02  respectively.  They  are  related  by  the 
pedaling  operation. 

in  another  higher  dimensional  surface  (f)2.  The  equation  of  motion  for  02  is  the  same 
as  the  standard  one  as  shown  in  (6.11),  but  we  do  not  need  to  solve  this  equation 
directly  for  function  02-  02  can  be  regarded  as  being  implicit  in  the  procedure,  we  do 
not  obtain  the  level-set  of  02,  namely,  the  snake  pedal  curve  Ve{t),  from  evaluating 
02  and  determining  the  points  02  =  0,  instead,  we  first  evaluate  0i  and  determine 
its  zero  set,  then  we  apply  the  pedaling  operation  to  obtain  the  snake  pedal  Ve{t). 


In  addition,  even  though  we  never  solve  02,  02  is  very  useful  in  the  derivation  of  the 
governing  PDE  for  the  function  0i.  The  development  of  the  relationship  between 
the  gradient  of  0i  and  02  is  a  crucial  step,  we  will  derive  this  relationship  later  in 
Section  6.4.3.  Note  that  the  "snake"  model  discussed  here  is  not  a  classical  energy- 
minimizing  model  any  more,  since  its  deformation  is  not  obtained  by  minimizing  some 
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of  clustered  obj 


internal  deformation  energy,  but  we  will  still  use  the  name  "snake"  to  represent  the 
controlling  active  contours  in  the  snake  pedal  model.  Fig.  6.2  depicts  this  important 
relationship  between  'P{t),Ve{t),(j)i,  and  (j)2. 

Note  that  an  important  feature  of  this  modeling  technique  is  the  incorpora- 
tion of  a  global  parameterized  shape,  namely,  the  generator,  into  the  curve  evolution 
framework.  This  global  parameterized  shape  can  be  very  useful  in  applications  in- 
volving recognition  of  clustered  objects.  For  example,  in  Fig.  (6.3),  there  are  three 
clusters  of  objects  in  2D.  Each  cluster  has  its  center,  size  and  orientation,  which  can 
be  represented  by  an  ellipse.  The  snake  pedal  model  can  capture  this  "global  orien- 
tation" of  the  cluster  via  the  global  part  of  the  model,  namely,  the  generator,  while 
recovering  every  single  object  in  the  cluster.  Since  a  geometric  active  contour  is  used 
to  represent  the  controlling  "snake"  in  the  model,  and  also  the  model  can  capture  the 
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"global  orientation"  of  a  group  of  objects,  this  snake  pedal  model  is  referred  to  as  a 
hybrid  geometric  active  contour  model.  Detailed  discussion  on  the  related  concepts 
will  be  given  in  the  following  sections,  starting  with  the  derivation  of  the  relationship 
between  the  snake  and  the  snake  pedal  evolution. 

6.4.2    Relation  Between  the  Snake  and  Snake  Pedal  Curve  Evolution 

Consider  a  snake  p  =  (pi,P2)^  and  a  generator  a  —  {ai,a2)^  with  normal 
=  (rii,n2)^,  where  p  and  a  are  related  by  the  radius-based  correspondence  as- 
sociating scheme  described  in  Chapter  3,  we  obtain  the  corresponding  snake  pedal 
Pe  =  {Pei,Pe2)^  via  the  following  pedaling  operation: 


here  we  use  the  modified  pedaling  operation  given  in  Chapter  3  for  the  purpose  of 
explanation.  Rewriting  the  Eq.  (6.17)  in  component  form,  we  obtain 


Pe 


(6.17) 


(6.18) 


or  simply. 


Pe  =  J4  P  -  b, 


(6.19) 


where 


3a 


1  + 


(6.20) 


7l{+  7I2 
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and 


b  = 


a 


1  •  ni  +  Qo  •  712 


(6.21) 


Taking  the  inverse  of  J^,  denoted  by  J^,  we  have 


—  2 


(6.22) 


We  should  remind  the  reader  that  by  using  the  association  scheme  described  in  Chap- 
ter 3,  the  forms  of  and  b  do  not  change  over  time,  but  their  contents  do  change. 
But,  for  the  following  reason,  we  can  still  treat  as  a  constant  matrix  and  b  as  a 
constant  vector  with  respect  to  time  when  evolving  p(t)  or  Pe(t).  First,  the  genera- 
tor does  not  evolve  over  time,  we  only  estimate  the  generator  parameters  iteratively; 
Second,  the  changes  of  and  b  are  very  small  compared  with  that  of  p  or  Pe. 
Therefore,  Eq.  (6.19)  can  be  explicitly  written  as  a  function  of  time: 


Taking  partial  derivative  on  both  sides  of  Eq.  (6.23)  with  respect  to  time  t  yields 


Pe(t)  =         P(^)  -  b. 


(6.23) 


dt 


dp{t) 


dt 


(6.24) 


or 


dpit) 


dPejt) 

dt 


(6.25) 


dt 


88 


Eq.  (6.24)  and  (6.25)  reveal  a  very  crucial  relationship  between  the  snake  and  snake 
pedal  evolution:  The  evolution  of  two  curves  are  linearly  related  by  a  Jacobian  matrix. 
We  remind  the  reader  that  the  relationship  between  the  snake  and  snake  pedal  is  not 
linear,  but  their  evolutions  over  time  are  related  by  a  linear  transformation!  This 
linear  transformation  between  the  evolution  of  the  snake  and  the  evolution  of  the 
snake  pedal  largely  simplifies  our  further  derivation  and  computation,  as  seen  in  the 
following  sections. 

6.4.3    PDE  for  the  Snake  and  Snake  Pedal  Evolution 

With  the  relation  between  the  snake  and  the  snake  pedal  curve  evolution  given 
by  (6.25),  we  now  derive  the  equation  of  motion  for  the  higher  dimensional  surface 
01,  in  which  the  snake  is  embedded.  As  discussed  in  Section  6.3.2,  the  zero  level  set 
of  the  higher  dimensional  surface  0i  is  represented  by  Eq.  (6.8), 


{V{t)  e  3f?'  : 


<f>{p,t)  =  0} 


(6.26) 


Differentiating  Eq.  (6.8)  with  respect  to  t  yields 


+  V(/)i  ~  =  0. 


(6.27) 


at 


Similarly,  for  the  higher  dimensional  function  (j)2,  we  have  the  following  formula 
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From  the  discussion  in  the  previous  section,  we  have 


d 


or 


(6.29) 


As  described  earher,  we  want  to  minimize  the  arc-length  of  the  snake  pedal,  which 
amounts  to  imposing  a  smoothness  constraint  on  the  snake  pedal  curve,  hence  we 
require  that 


dt 


(6.30) 


where  F{ke)  is  the  speed  function,  which  depends  on  the  curvature  of  the  snake  pedal 
kg.  Ne  is  the  unit  normal  of  the  snake  pedal.  Similar  to  the  relation  between  the 
normal  of  the  snake  and  its  higher  dimensional  function  (/),  :  N  =  -m^^,  for  the 

iiv0iir 

zero  level  set  of  02,  the  following  relation  holds: 


IIV02 


(6.31) 


Combining  Eq.  (6.28),  (6.30)  and  (6.31),  we  obtain 


d(t>2 
dt 


=  F{k,)\\Vci>2\ 


(6.32) 
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Eq.  (6.32)  is  the  standard  level-set  evolution.  For  the  governing  equation  of  the  snake 
curve,  we  substitute  Eq.  (6.29),  (6.31),  and  (6.30)  into  Eq.  (6.27)  to  obtain 


^  +  V01  •  {JBF{ke)Ne)  =  0  (6.33) 


Substituting  Ne  in  Eq.  (6.33)  with  Eq.  (6.31)  yields  the  following  expression 
For  simplicity,  Eq.  (6.34)  can  be  rewritten  as 


Eq.  (6.32)  and  (6.35)  are  the  equations  of  motion  for  0i  and  (^2,  respectively.  We 
rewrite  them  together  as  follows: 


(6.36) 


Eq.  (6.36)  represents  a  coupled  equation  for  the  evolution  of  </>i  and  (f)2.  One  approach 
to  solve  the  PDE  system  (6.36)  is  to  use  a  combination  of  central  differences  and  the 
finite  upwind  diflFerence  method  [52]  to  solve  the  first  equation  in  every  iteration,  and 
then  solve  the  second  equation  using  a  similar  method  subsequently.  We  will  discuss 
the  up-wind  finite  difference  method  in  Chapter  7.  This  approach  is  straightforward, 
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but  needs  the  memory  storage  for  both  0i  and  02  at  every  grid  point.  Note  that 
in  Eq.  (6.36),  only  the  gradients  of  01  and  02  are  involved  on  the  right  hand  sides 
of  both  equations,  we  propose  a  more  elegant  approach  to  solve  the  PDE  system  of 
(6.36)  simultaneously  with  much  less  storage  requirement,  by  employing  the  intrinsic 
relation  between  V0i  and  V02  as  discussed  in  the  next  section. 


6.4.4    Relation  Between  V0i  and  Vc!>9 

To  discover  the  relation  between  V0i  and  V02,  at  first,  we  arbitrarily  choose 
representative  points  on  the  snake.  Let  p  =  (x,  yY  denote  the  coordinates  of  point 
a.  The  other  two  representative  points  h  and  c  are  chosen  by  shifting  point  a  in  x 
direction  and  in  y  direction,  such  that  p  =  {x+/S.x,  yY  and  p  =  (x,  y+/S.yY  are  their 
coordinates,  respectively.  Let  a', 6',  and  c'  be  the  corresponding  points  on  the  snake 
pedal  that  are  obtained  by  applying  the  pedaling  operation  with  respect  to  point 

a,6,  and  C,  with  coordinates  p^  =  {Pex^PeyV ,  Pe  =  {Pex,PeyY  and  p^  =  {Pex,PeyV^ 

respectively.  More  specifically, 


(6.37) 


Pex  — 


Pex{x  +  Ax,y) 


< 


(6.38) 


Pey  — 


Pey{x  +  Ax,y), 


and 


(6.39) 
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(x,y+Ay) 


(x,y) 


(x+Ax,y) 


a,  b,  c  :    points  on  the  snake 

a',  b',  c' :  points  on  the  snake  pedal 


(a) 


Figure  6.4.  Three  representative  point-pairs  on  the  snake  and  snake  pedal. 


The  higher  dimensional  function  values  at  points  a,  b,  and  c  are  given  by  0i,  4>i 
and  4>i  respectively.  We  construct  the  higher  dimensional  surfaces  (/>i  and  ^2  in  such 
a  way  that  the  zero  level  sets  of  (pi  and  02  represent  the  snake  and  the  corresponding 
snake  pedal  curve,  respectively.  The  same  relation  holds  for  the  other  level  sets  of 
(^1  and  (j)2.  Therefore,  the  higher  dimensional  function  values  of  (f)2  at  position  a',  b' 
and  c'  are  0i,  ^1,  and  4>i,  respectively. 

Since  (^,  ^)'^  is  the  gradient  of  </)2  at  location  a',  from  location  a'  to  b', 
using  the  definition  of  directional  derivatives,  we  have  the  first  expression  in  Eq. 


93 


(6.40)  and  from  a'  to  c'  the  second  relation  holds: 


df2 
dx 


(Pex-Pex)     +     ^{Pey-Pey)     =    </^l  "  01 
Pex)     +     ^(Pey-Pey)     =  ^1-01- 


(6.40) 


Similarly,  given  the  gradient  of  <pi,  namely,  (^^,  ^^^^  location  a,  from  a  to  b, 
we  have 


(6.41) 


Since  the  snake  and  the  snake  pedal  are  related  by  the  pedaling  operation  in 
Eq.  (6.39),  using  the  differentiation  theory  of  function  of  several  variables,  we  have 


Pex-Pex     =  ^Ax+^Ay 


^Pey^^  ,  dPsJL 


Pey      Pey     —  Ax 


Ay, 


(6.42) 


Since  from  point  a'  to  6',  Ay  =  0,  Eq.  (6.42)  can  be  simplified  as: 


Pex      Pex     —  Ax 


dp, 


.  Pey      Pey  - 


^Ax. 


(6.43) 


Similarly,  from  a'  to  c',  we  have: 


Pex  -  Pex 
Pey  ~  Pey 


^Ax  +  ^Ay 
^Ax  +  ^Ay, 


(6.44) 
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Since  Ax  =  0  from  a'  to  d  , 


Pex-Pex  = 


Pey      Pey  - 


Ay. 


(6.45) 


Substituting  (6.41),  (6.43),  and  (6.45)  in  to  (6.40)  yields: 

'  ^  .  dp^Ax  +  ^  •  ^Ax  = 
I    ox     ox  oy  ox 


OX 


dd 


dy 


(6.46) 


Eliminating  Ax  and  Ay  on  both  sides  of  Eq.  (6.46)  results  in 


dy     dx  dx 

ex       ,       902  .  dPey      ^  90^ 

dy^  ~By  'By 


ox  ox 

d^.dp, 
ox  oy 


(6.47) 


Writing  in  the  matrix  form,  we  obtain 


901 

dpex 

dx 

OX 

d(l>i 

dpex 

L  dy  J 

L  dy 

dPey 

dx 

902 

dx 

dpey 

dy  J 

902 

L  dy  J 

(6.48) 


or  simply, 


V01  =  3a  ■  V02, 


(6.49) 


or 


V02  =  3  b  ■  y(f>u 


(6.50) 
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where  the  Jacobian  between  the  evolution  of  the  snake  and  snake  pedal  curve  and 
Jfi  =  J^^  are  defined  in  Eq.  (6.20)  and  Eq.  (6.22),  respectively. 

After  a  simple  yet  tedious  derivation,  we  obtain  a  surprisingly  simple  relation 
between  the  V0i  and  V(?!>2.  At  first  glance,  this  relation  is  seemingly  contradictory  to 


between  the  gradients  of  0i  and  (j)2,  the  relation  should  be  V02  =  J>iV^i,  rather  than 

*    t  ■  * 

The  development  of  the  relationship  between  V(/»i  and  V02  is  a  crucial  step  in 
the  derivation  of  the  evolution  of  (pi ,  which  will  be  given  in  the  next  section. 

6.4.5    PDE  for  the  Higher  Dimensional  Function 

Since  our  objective  is  to  obtain  the  equation  of  motion  for  the  higher  dimensional 
surface  (pi,  we  substitute  (f)2  with  Jb  0i  in  Eq.  (6.35)  to  get 


Note  that  we  use  the  symmetric  property  of  Jb  in  the  above  derivation  from  the 
second  step  to  the  third  step.  The  representation  of  kg  in  terms  of  (pi  is  quite  com- 
plicated, we  refer  the  reader  to  Appendix  B  for  these  details. 

Eq.  (6.51)  is  the  equation  of  motion  for  (pi  and  constitutes  the  primary  equation 
in  our  application.  When  using  the  snake  pedal  model  for  recovery  of  the  shape  of 


our  imagination.  One  would  think  that  since 


5^  =  J^^l^,  if  there  is  any  Jacobian 


(6.51) 
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interest  from  image  data,  we  need  to  solve  the  more  general  equation  of  motion  for 
the  higher  dimensional  function  (^i.  From  Eq.  (6.15),  by  using  a  similar  approach, 
we  obtain  the  more  general  equation  of  motion: 

^  =       y)  F{K)  II  Jb V<^i  II  +  Vc^i  •  Vp,  (6.52) 

where  g{x,  y)  is  an  image  feature  based  function  and  is  used  to  stop  the  curve  evolution 
when  the  contour  is  close  to  the  desired  edges.  V</»i  and  must  be  evaluated  at 
locations  on  the  snake  and  the  snake  pedal,  respectively. 


CHAPTER  7 

NUMERICAL  SOLUTIONS  FOR  SHAPE  RECOVERY 
WITH  HYBRID  GEOMETRIC  ACTIVE  CONTOURS 


In  Chapter  6,  we  discussed  that  how  the  snake  pedal  may  be  used  to  recover 
the  object  boundaries  from  an  image  using  a  novel  global  plus  local  shape  estimation 
procedure.  For  the  global  shape  estimation,  we  employ  the  Levenberg-Marquardt 
(LM)  method;  while  for  the  local  shape  estimation,  we  adopt  a  modified  level-set 
method.  We  have  discussed  the  LM  method  in  Chapter  5,  we  will  now  discuss  the 
modified  level-set  framework  for  local  shape  recovery. 

In  Chapter  6,  we  illustrated  that  the  solution  for  the  snake  pedal  evolution 
can  be  achieved  by  first  solving  the  snake  evolution,  which  is  embedded  in  a  higher 
dimensional  surface  (^i,  then  applying  the  pedaling  operation  on  the  snake.  There- 
fore, developing  reliable  and  efficient  numerical  algorithms  for  solving  the  governing 
equation  of  that  is,  Eq.  (6.51),  becomes  one  of  the  main  tasks  in  applying  the 
snake  pedal  model  for  extracting  shapes  of  interest  from  image  data.  The  governing 
equation  of  motion  for  0i,  in  the  simplest  form  is  given  in  Eq.  (6.51),  we  rewrite  this 
in  the  following  for  the  purpose  of  discussion, 

^  =  F(A;e)  ||JbV</.x||,  ,  (7.1) 
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where  Js  is  the  Jacobian  between  the  snake  and  snake  pedal  evolution.  In  our 
implementation,  we  assume  that  is  a  constant  matrix  with  respect  to  time  based 
on  the  reasons  stated  in  Section  6.4.2.  kg  is  the  curvature  of  the  snake  pedal,  and 
F{.)  is  the  speed  function.  To  discuss  the  numerical  behavior  of  (7.1),  we  will  first 
discuss  the  solution  for 


Eq.  (7.2)  is  equivalent  to  (7.1)  if  Jb  is  an  identity  matrix  and  kg  is  equal  to  k,  the 
curvature  of  the  snake.  Eq.  (7.2)  is  the  standard  level-set  equation  and  its  solution 
has  been  widely  discussed  in  literature  [52,  62].  We  will  show  that  the  solution 
technique  for  (7.2)  with  minor  modifications  can  be  apphed  to  solve  (7.1).  We  will 
then  present  some  shape  recovery  results  by  applying  the  snake  pedal  model  to  image 
data  in  Chapter  8. 

7.1    Numerical  Algorithm  for  Standard  PDE-based  Curve  Evolution 

In  Chapter  6,  we  discussed  the  eff'ect  of  the  advection  term  Fa  and  the  diffusion 
term  Fg  in  the  speed  function  F{k).  Fig.  7.1  depicts  the  curve  evolution  with  two 
different  speed  functions.  In  Fig.  7.1  (a),  the  curve  moves  with  unit  speed,  while 
in  Fig.  7.1  (b),  the  speed  depends  on  the  curvature  of  the  curve.  To  illustrate  the 
numerical  behavior  of  (7.2),  we  substitute  F{k)  =  1  -  eA;  as  a  typical  speed  function 
and  the  equation  of  motion  becomes 


F{k)  \\Vcf> 


(7.2) 


dt 


{Fa  +  Fg)  \\\/cf>,\\  =  {l-ek) 


dt 


(7.3) 
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Constant  speed  Curvature  dependence 

F(k)=l  F(k)=l-  £k 

(a)  (b) 


Figure  7.1.  Curve  Evolution:  The  curve  moves  with  unit  speed  in  (a),  and  moves 
with  curvature  dependent  speed  in  (b).  The  soUd  Unes  are  the  initial  curves  and  the 
dashed  lines  are  the  curves  after  evolution. 

where  we  use  Fa  =  1  and  Fq  =  -ek.  Eq.(7.3)  poses  an  initial  value  problem: 

0u  +  (<^L  +  0U^/^  =  eV-(^),  (7.4) 
with  the  initial  condition: 

(l)i{^,y,i^O)^±d{x,y),  (7.5) 

where  d{x,  y)  is  the  distance  from  {x,  y)  to  the  initial  snake  curve  'P(O).  As  discussed 
in  Sethian  [62],  with  the  presence  of  the  diffusion  term  Fq  (curvature  dependent 
term),  for  e  >  0,  the  parabolic  right-hand  side  diffuses  sharp  gradients  and  force 
to  stay  smooth  for  all  values  of  t.  When  there  is  no  diffusion  term  F^,  e  =  0,  the 
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curve  moves  with  unit  speed,  and  a  corner  may  develop  from  smooth  initial  data. 
Once  a  corner  develops,  it  is  not  clear  how  to  propagate  the  curve  in  the  normal  di- 
rection, since  the  derivative  is  not  defined  there.  A  variety  of  "weak"  solutions  have 
been  proposed  to  evolve  the  curve  beyond  the  formation  of  singularities.  Of  all  such 
weak  solutions,  one  is  interested  in  the  one  that  is  the  limit  of  smooth  solutions  as 
e  ->  0.  Another  way  to  characterize  this  weak  solution  is  via  the  following  "entropy 
condition"  [52]:  if  the  curve  is  viewed  as  the  boundary  of  a  propagating  flame,  then 
once  a  particle  is  burnt,  it  stays  burnt.  Hence,  approximations  to  the  spatial  deriva- 
tives are  sought  that  do  not  artificially  smooth  sharp  corners  but  pick  out  the  correct 
entropy  solution  beyond  the  occurrence  of  singularities.  Careful  adherence  to  this 
stipulation  produces  the  Huygens  principle  construction.  In  fact,  the  schemes  given 
by  Osher  and  Sethian  [52,  62]  are  motivated  by  the  fact  that  the  entropy  condition 
for  evolving  curves  is  identical  to  the  one  for  hyperbolic  conservation  laws,  where 
stable,  consistent,  entropy-satisfying  algorithms  have  a  rich  history  [22,  33,  40,  41]. 

Let's  consider  the  curve  evolution  in  one  dimensional  space  with  constant  normal 
velocity  F  =      =  1,  e  =  0.  Eq.  (7.1)  becomes  a  standard  Hamilton- Jacobi  equation 

0u  +  i^(0ix)  =  0,  (7.6) 

with  H{(j)ix)  =  -{(klxY'"^  and  a  given  initial  value  of  (f)^.  Let  u  -  </)i^,  by  taking  the 
derivative  with  respect  to     Eq.  (7.6)  becomes 


+  [E{u%  =  0, 


(7.7) 
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with  H{u)  —  -{u^yi'^.  Eq.  (7.7)  is  a  scalar  hyperbolic  conservation  law  in  one  space 
variable.  Numerous  research  results  reported  in  literature  [52,  61]  have  demonstrated 
that  if  we  use  the  standard  forward  difference,  central  difference,  or  backward  differ- 
ence mechanisms  to  discretize  this  equation,  the  solutions  can  develop  discontinuous 
jumps,  known  as  shocks,  even  with  smooth  initial  data. 

As  in  Osher  and  Sethian  [52],  we  use  a  conservative,  monotone  scheme  called 
upvirind  finite  diff"erence  method  to  discretize  Eq.(7.7).  The  salient  point  to  be 
noted  is  that  this  conservative  scheme  produces  a  solution  that  satisfies  the  entropy 
condition.  An  upwind  scheme  automatically  differences  in  the  outward-flowing  di- 
rection during  the  curve  evolution,  thus  naturally  accommodates  for  the  evolution 
behavior.  We  now  introduce  the  upwind  finite  difference  scheme  to  solve  Eq.  (7.6). 

Writing  Eq.  (7.6)  using  a  forward  difference  in  time  yields: 

<i>lt'  =  K-^tH{ul  (7.8) 

where  (l)^^  denotes  the  value  of  0i  at  grid  location  i  Ax  and  time  n  At,  with  Ax 
and  At  being  space  and  time  intervals  in  the  finite  difference.  In  the  upwind  finite 
difference,  we  use  an  appropriate  numerical  flux  function  G  to  approximate  H 

rit'     =  AiG(u,_i/2,n,+i/2) 

(7.9) 

=  4>l-^tG{D-cl>,,,Dtci>,,), 
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where       and       are  standard  forward  and  backward  difference  operators: 

C;^,  =  «^t2  (7.10) 

Dth.  =  (7.11) 

We  can  use  the  Hamilton-Jacobi  fiux  function  given  in  Osher  and  Sethian  [52]  as  the 
flux  function  G.  When  H{u)  can  be  written  as  a  function  of  u^,  that  is,  H{u)  —  f{u^) 
for  some  function  /,  the  flux  function  is: 

G(Mi_i/2,u,+i/2)   =  G{D-(j)u,D+(f)u) 

(7.12) 

=   f{{max{D-cl>u,  0)f  +  {min{D+<j>u,  O))^). 

This  conservative  monotone  scheme  is  an  upwind  method  in  that  it  differences  in  the 
direction  of  propagating  characteristics.  In  the  case  when  Fa  =  1,  /(u^)  =  -(w^)^/^, 
Eq.  (7.9)  reduces  to 

=  <I>1  -  At{  {max{D-^u,  0))'  +  {min{Dthi,  0))'  Y'\  (7.13) 

This  algorithm  produces  the  correct  entropy-satisfying  weak  solution  to  the  curve 
evolution  problem.  ' 

7.2    Our  Numerical  Approach 

The  above  discussion  can  be  extended  to  problems  in  higher  dimensions  [45,  52]. 
Recall  that  the  original  intent  was  to  solve  (7.3).  In  two  dimensions,  the  scheme  given 
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in  Eq.  (7.13)  is  extended  by  differencing  in  each  direction  to  produce  the  following 
scheme  for  Eq.  (7.3) 

x{(maa;(D-0i,(ij),O))2  +  (mm(D+</)i,(ij),  0))^+ 

(7.14) 

(max(/);  0))2  +  (m2n(D+</.i,(,,),0))2}^/^ 

-AtFc  ||V0i||, 

where  <Ai,(ij)  represents  the  value  of  at  grid  location  (zAx,jA?/)  and  time  nAt. 
Ax  and  Ay  are  space  intervals  in  the  finite  differences  along  the  x,y  directions,  and 
At  is  the  time  interval.  Since  the  curvature  dependent  term  is  a  parabolic  term, 
standard  central  difference  approximation  can  produce  stable  solutions  for  this  term, 
no  special  treatment  needs  to  be  done  for  Fg||0i||  =  The  same  conclusion 

holds  if  we  substitute  khy  ke,  the  curvature  of  snake  pedal. 

To  solve  the  original  equation  of  motion  for  the  snake:  =  F{ke)  ||JbV(/)i||, 
for  the  diffusion  term,  we  use  the  standard  central  difference  approximation;  For  the 
advection  term,  we  need  to  solve  the  following  hyperbolic  initial-value  problem: 

(t>u  =  \J a{x,  y)(f>l^  +  b{x,  y)(j)'iy  +  c{x,  y)(f)ix(t>iy,  (7.15) 

where  a{x,y),b{x,y)  and  c{x,y)  are  determined  by  the  entries  of  3b  and  do  not 
change  over  time.  and  (f^jy  can  still  be  approximated  by  the  upwind  finite  dif- 
ference scheme  discussed  earlier.  But  for  the  <f)u(l>iy  term,  we  use  the  minmod 
finite  difference  approximation  discussed  by  Kimmel  et  al.    [34].   The  minmod 
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finite  derivative  is  defined  as: 


sign{a)min{\a\,  \b\)   if  a6  >  0 
minmod{a,b}  =  I  (7-16) 

0  otherwise. 


We  use  this  definition  to  approximate  (j)ix4>iy  by 

<t>\x(t>\y\x=iAx,y=i^y  =  minmod{D^ (pi^^ij),  D' (j)^^^,^^)  minmod{D^(j)^^(i^j),Dy(j)x,{i,j)), 

(7.17) 

where  D^,D~,D^,  and  D~  are  as  defined  before. 

Combining  these  finite  differences  yields  a  first  order  numerical  scheme  for  solv- 
ing (7.15),  the  advection  term  in  Eq.  (7.1): 

[aij((maa;(L»-(^i,(,j),0))2  +  (mzn(L)+0i,(,j),  0))^) 
+6ij((maa;(D;0i,(i,,),  0))^  +  (mm(L>+</.i,(.,,),  0))^) 
-Ci,jmmmorf(D+0i,(ij),D-0i,(ij))mmmocfe(D+^ 

(7.18) 

This  numerical  scheme  is  stable  and  provides  a  natural  way  to  handle  the  topological 
changes  in  the  snake  evolution.  For  numerical  schemes  with  higher  order  of  accuracy 
that  deal  with  such  Hamilton-Jacobi  type  of  equations,  we  refer  the  reader  to  Osher 
and  Shu[68]. 


CHAPTER  8 

SHAPE  RECOVERY  EXPERIMENTS  WITH  SNAKE  PEDALS 


This  chapter  presents  2D  and  3D  shape  recovery  examples  using  the  snake  pedal 
model  developed  in  this  thesis.  Examples  include  shape  recovery  from  a  variety  of 
2D/3D  data  sets,  including  range  data,  MRI  and  CT  image  data.  The  examples 
shown  depict  a  wide  range  of  shapes  with  complex  geometries  and  topologies.  These 
experiments  serve  to  demonstrate  the  power  and  versatility  of  the  snake  pedal  model 
for  shape  recovery  and  shape  synthesis. 

8.1    Shape  Recovery  with  Snake  Pedals  in  2D 

In  this  section,  we  present  a  set  of  three  shape  recovery  experiments  in  2D. 
All  the  experiments  are  performed  on  MRI  brain  scans.  In  all  the  experiments,  the 
spatial  resolution  of  the  images  is  256  x  256  and  a  region  of  interest  is  identified 
interactively. 

In  the  first  experiment,  we  show  shape  estimation  of  a  caudate  nucleus,  a  cor- 
tical structure  in  the  human  brain,  from  a  slice  of  a  brain  MRI.  The  snake  pedal 
initialization  is  shown  in  Fig.  8.1  (a),  while  (b)  and  (c)  show  the  intermediate  stage 
and  final  model  fit  respectively.  This  fitting  is  achieved  at  interactive  rates  on  an 
SGI-02.  Fig.  8.2  depicts  a  similar  experiment  carried  out  on  a  cerebellum  structure 
in  the  human  brain,  from  a  slice  of  brain  MRI.  Fig.  8.2  (a)  -  (c)  show  the  snake 
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(a)  (b)  (c) 

Figure  8.1.  Fitting  result  (c)  of  a  snake  pedal  to  a  caudate  nucleus  with  initialization 
shown  in  (a)  and  intermediate  stage  shown  in  (b). 

pedal  initialization,  intermediate  stage,  and  final  fits  respectively.  Fig.  8.3  (a)  -  (c) 
depict  a  2D  fitting  example  for  a  hippocampus  in  the  human  brain.  Fig.  8.3  is  the 
initialization,  (b)  and  (c)  are  the  intermediate  stage  and  the  final  fit  respectively. 
In  the  caudate  nucleus  and  the  cerebellum  examples,  we  use  both  image  gradient 
potential  as  well  as  spring-based  potential  from  the  data  points  as  the  external  force 
to  guide  the  model  to  conform  to  the  desired  image  features.  In  the  hippocampus 
example,  we  only  use  the  spring-based  data  point  force,  since  the  image  boundary  is 
not  well  defined.  The  points  (shown  in  red)  are  placed  by  an  expert  neuroscientist. 

Note  that  the  model  initialization  is  quite  far  from  the  final  true  shape  in  each 
example.  Also,  the  model  fits  are  quite  accurate  qualitatively  (for  visual  inspection). 
The  accurate  extraction  of  the  structures  in  the  human  brain  plays  an  important  role 
for  the  computer  aided  medical  diagnosis.  For  example,  when  a  person  is  afflicted  with 
epilepsy,  the  volume  of  the  hippocampus  in  the  person's  brain  changes  significantly 
from  the  volume  of  the  hippocampus  in  the  normal  brain  [21,  39]. 
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(a)  (b)  (c) 

Figure  8.3.  Fitting  result  (c)  of  a  snake  pedal  (in  green)  to  a  hippocampus  with 
initialization  shown  in  (a)  and  intermediate  stage  shown  in  (b).  The  data  points  (in 
red)  are  placed  by  the  expert  along  the  hippocampus  boundary. 
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8.2    Shape  Recovery  with  Snake  Pedals  in  3D 

In  this  section,  we  present  a  set  of  seven  experiments,  three  of  these  experiments 
demonstrate  the  shape  reconstruction  with  the  snake  pedal  model  for  medical  image 
analysis  applications.  The  third  and  fourth  examples  are  designed  to  illustrate  the 
global  deformations  of  the  model  without  explicitly  including  the  bending,  tapering, 
and  twisting  parameters  in  the  model  representation.  In  these  two  examples,  the 
model  is  fitted  to  3D  synthetic  data  from  a  bent  shape  (hair  pin)  and  a  twisted 
shape,  respectively.  The  sixth  and  seventh  examples  demonstrate  model  fitting  to 
3D  sparse  range  data  points. 

The  medical  examples  are  organized  as  follows:  in  Fig.  8.4,  left  to  right,  the  im- 
ages depict  a  slice  of  an  MR  brain  scan  in  which  the  shape  of  interest — hippocampus, 
a  gyrus,  and  a  cerebellum — has  been  identified  by  a  neuroscientist  via  sparsely  placed 
points  (only  in  selected  slices  of  the  MR  scan)  on  the  shape  boundary.  The  next  image 
shows  the  collection  of  these  3D  points  (in  red)  and  the  initialized  snake  pedal  model 
(in  green)  followed  by  images  depicting  an  intermediate  stage  of  fitting  and  the  final 
fitted  model  respectively.  As  evident,  the  model  achieves  a  visually  accurate  fit  to 
the  data.  In  the  case  of  a  hippocampus,  model  fitting  was  performed  for  30  left  and 
right  brains  and  computed  volumes  were  validated  against  manual  measurements. 
The  Table  8.1  depicts  the  ratio  of  the  left  to  right  hippocampal  volumes.  The  degree 
of  match  between  the  manual  and  computed  volume  ratios  is  very  encouraging  and 
sets  the  stage  for  further  research  in  automatic  methods  for  epilepsy  studies. 
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(a) 

(b) 

(c) 

(d) 

(e) 

(f) 

(g) 

(h) 

• 

(i) 

0) 

(k) 

(1) 

Figure  8.4.  Fitting  results  of  snake  pedal  surfaces  to  a  hippocampus  (d),  a  gyrus  (h), 
and  a  cerebellum  (1)  with  initializations  shown  in  (b),  (f),  and  (j),  the  intermediate 
stages  are  shown  in  (c),  (g),  and  (k),  respectively,  (a),  (e)  and  (i)  are  slices  of  MRI 
scans  for  the  hippocampus,  gyrus,  and  cerebellum,  respectively. 
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Table  8.1:  Comparison  between  the  manual  and  computed  ratios 


of  the  left  to  right  hippocampal  volumes 


Manual 

Automatic 

Index 

DilllU 

Volume 

L/R 

Volume 

L/R 

4897  R 

19 

3.3345 

3.4142 

1 

4897  L 

lb 

3.3310 

0.9895 

3.6625 

0.9322 

4915  R 

20 

3.8114 

3.7209 

2 

4915  L 

19 

3.7668 

0.9882 

3.6412 

0.9785 

4957  R 

19 

3.8114 

3.7209 

3 

4957  L 

20 

3.7668 

0.9882 

3.6412 

0.9785 

4816  R 

11 

3.5521 

3.5430 

4 

4815  L 

OA 
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The  two  synthetic  data  examples  are  organized  as  follows:  in  Fig.  8.5,  left 
to  right,  the  images  depict  a  3D  synthetic  data  set  (in  red)  along  with  the  model 
initialization  (in  green)  superimposed,  which  is  followed  by  images  of  intermediate 
stage  of  fitting  and  the  final  model  fit.  Note  that  in  these  synthetic  data  examples, 
in  one  case,  the  data  contains  shapes  exhibiting  large  bending  and  in  another  a  large 
amount  of  twist.  What  we  would  like  to  stress  here,  is  the  fact  that  the  model  did 
not  require  any  explicit  bending  or  twisting  transformations  (unlike  the  deformable 

t  4  •. 
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(d)  (e)  (f) 


Figure  8.5.  Fitting  results  of  snake  pedal  surfaces  to  a  bend  (c)  and  twist  (f)  shape. 
Initializations  are  shown  in  (a)  and  (d)  followed  by  the  intermediate  stages  (b)  and 
(e)  respectively. 

superquadric)  to  achieve  the  fits.  As  mentioned  earlier,  modeling  schemes  that  require 
explicit  incorporation  of  bending,  twisting  transforms  are  highly  nonlinear  and  are 
known  to  introduce  numerical  instabilities  into  the  fitting  algorithms.  In  contrast,  our 
modeling  scheme  does  not  require  explicit  application  of  such  bending  and  twisting 
transformations  to  achieve  the  desired  fit. 

The  sixth  and  seventh  examples  depict  the  snake  pedal  fitting  in  3D  to  a  "vul- 
cau"  (Spock  form  Star  Wars)  head  and  a  mechanical  part  (anvil)  from  range  data. 
The  images  in  Fig.  8.6  are  arranged  in  the  same  way  as  the  examples  in  Fig.  8.2. 

In  all  the  examples,  we  use  a  20  x  40  mesh  to  discretize  the  model  except  in  the 
bent  shape  example  (Fig.  8.5),  in  which  we  use  a  40  x  80  mesh. 
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Figure  8.6.  Fitting  results  of  snake  pedal  model  to  spock(c)  and  an  anvil  (f).  Initial- 
izations are  shown  in  (a)  and  (d)  followed  by  the  intermediate  stages  of  fitting  in  (b) 
and  (e)  respectively. 
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(a)  (b)  (c)  (d) 


Figure  8.7.  Hybrid  geometric  active  contour  model  fit  examples:  (a)  is  model  ini- 
tializations (snake  pedal  in  green  and  generator  in  red),  (b)  and  (c)  are  intermediate 
stages  of  evolution  and  final  fits  are  shown  in  (d). 

8.3    Fitting  Results  with  the  Hybrid  Geometric  Active  Model 

In  this  section,  we  present  the  level-set  implementation  of  the  hybrid  geometric 
active  model  with  an  ellipse  generator  described  in  Chapter  6.  In  Fig.  8.7  and 
Fig.  8.8,  left  to  right,  images  show  model  initialization — depicting  the  snake  pedal  in 
green  and  the  global  shape  in  red,  intermediate  stages,  and  final  fits,  respectively.  The 
human  head  image  is  obtained  from  a  slice  of  an  MRI  scan,  and  the  spanner  image 
is  obtained  from  the  "Computer  Vision  Test  Images"  web  site  [1].  As  is  evident,  the 
snake  pedal  was  successful  in  extracting  the  shape  of  interest  in  these  images  via  the 
model  fitting  process  wherein  an  image-based  speed  function  was  used  to  deter  the 
model  evolution.  Note  that  the  "global"  shape  depicted  in  these  examples  are  scaled 
versions  of  the  generator.  The  "global"  shape  provides  us  with  a  "global  size"  and 
"global  orientation"  of  the  shape  of  interest. 

More  medical  data  fitting  examples  are  given  in  Fig.  8.9.  In  each  row  of  the 
figure,  from  left  to  right,  images  show  the  model  initializations,  intermediate  stages 
and  the  final  model  fits.  The  first  and  second  examples  contain  the  2D  slices  from 


(a)  (b)  (c)  (d) 

Figure  8.8.  Hybrid  geometric  active  contour  model  fit  examples:  (a)  is  model  ini- 
tializations (snake  pedal  in  green  and  generator  in  red),  (b)  and  (c)  are  intermediate 
stages  of  evolution  and  final  fits  are  shown  in  (d). 

an  abdominal  CT  scan.  In  the  first  row,  the  model  is  initialized  in  the  stomach  and 
in  the  second  row,  the  model  is  initialized  inside  the  liver  metastasis.  The  third  and 
fourth  example  contain  2D  slices  from  an  MR  scan  of  the  human  heart.  The  model 
is  initialized  to  capture  the  endocardium  and  the  myocardium  structure  respectively. 
The  snake  pedal  is  shown  in  green  and  the  global  shape  in  red.  We  use  image-based 
speed  function  to  deter  the  model  evolution  in  all  the  four  examples. 

Figure  8.10  depicts  four  topological  change  examples.  In  each  row,  from  left  to 
right,  images  depict  model  initialization,  intermediate  stages  of  evolution  and  final  fit 
respectively.  In  the  first  example,  the  concept  of  a  "global"  shape  may  be  associated 
with  a  "global"  orientation  of  a  cluster  of  shapes  (in  this  case,  two).  In  the  second 
and  fourth  example,  the  snake  pedals  are  initialized  to  include  all  the  objects  in  the 
images,  then  the  models  shrink  and  split  to  fit  to  all  the  object  boundaries.  While  in 
the  third  example,  the  snake  pedal  is  initialized  as  a  single  small  ellipse,  as  the  fitting 
proceeds,  the  model  expands  and  splits,  and  finally  fits  to  all  the  object  contours  in 
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Figure  8.9.  Hybrid  geometric  active  contour  model  fit  examples:  Left  to  right,  model 
initialization  (snake  pedal  in  green  and  generator  in  red),  intermediate  stages  of 
evolution  and  final  fit. 


118 

the  whole  image.  The  global  shapes  of  the  generator  are  not  shown  here  since  the 
meaning  of  the  "global"  shape  in  these  examples  is  not  very  useful. 

In  all  these  examples,  we  implemented  the  level-set  form  of  the  equation:  = 
g{x,y)  F{k,)  ||JbV</.i||  +  V<^i  ■  Vg,  Where,  g{WI)  =  1/(1  +  \\V{G  *  I)\\yK)  with 
G  *  I  being  a  Gaussian  convolved  with  the  image  and  K  being  a  scaling  constant. 
We  use  the  upwind  difference  and  minmod  difference  method  described  in  Chapter  7 
to  implement  this  Hamilton-Jacobi  equation  of  motion. 
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Figure  8.10.  Topological  change  examples  with  hybrid  geometric  active  contour  mod- 
els: Left  to  right,  model  initialization  (snake  pedal  in  yellow  or  green  and  generator 
in  red),  intermediate  stages  of  evolution  and  final  fit. 


CHAPTER  9 
CONCLUSIONS  AND  FUTURE  DIRECTIONS 


In  this  thesis,  we  have  presented  a  powerful  geometric  shape  modeling  scheme 
which  allows  for  the  representation  of  global  shapes  with  local  detail  and  permits 
model  shaping  as  well  as  topological  changes  via  physics-based  control.  Our  mod- 
eling scheme  consists  of  representing  shapes  by  pedal  curves  and  surfaces.  Physics- 
based  control  was  introduced  for  shaping  these  geometric  models  by  using  a  dynamic 
spline  to  represent  the  position  of  the  pedal  point  in  the  regular  pedal  curve/surface. 
The  model,  dubbed  as  a  "snake  pedal,"  allows  for  interactive  manipulation  by  ap- 
plying forces  to  the  snake.  A  level-set  implementation  of  the  geometric  active  con- 
tour facilitates  topological  changes  in  the  snake  pedal  model.  The  hybrid  geometric 
active  contour  model  in  essence  introduces  a  vehicle  for  incorporating  a  global  pa- 
rameterized shape  descriptor  into  the  now  popular  framework  of  geometric  active 
contours  /  surfaces . 

9.1    Salient  Features  of  the  Model 

The  snake  pedal  models  are  bestowed  with  several  unique  features  which  make 
them  very  powerful  and  attractive  for  use  in  shape  synthesis  and  recovery  applica- 
tions. A  list  of  these  features  are: 
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9.1.1  Geometric  Model  with  Physics-based  Control 

The  deformable  superquadric  model  is  created  by  superimposing  a  dynamic 
spline  on  a  geometric  primitive,  namely,  a  conventional  superellipsoid.  The  snake 
pedal  model  is  created  by  applying  the  pedaling  operation  to  a  geometric  primitive 
with  respect  to  a  dynamic  spline.  Both  models  allow  for  the  representation  of  global 
and  local  shape  characteristics  of  an  object,  but  the  snake  pedal  model  is  inherently  a 
geometric  model,  while  the  deformable  superquadric  model  is  a  physics-based  model. 

9.1.2  Compact  Representation 

Unlike  the  deformable  superquadric  model,  the  snake  pedal  model  does  not  need 
to  explicitly  include  the  bending,  tapering  and  twisting  parameters  in  its  represen- 
tation to  realize  the  global  deformation  for  the  same.  Since  there  are  no  additional 
nonlinearities  introduced  into  the  model,  numerical  efficiency  and  stability  can  be 
largely  improved  in  the  model  fitting  to  the  image  data  with  the  snake  pedal  model. 

9.1.3  Automatic  Handling  of  Topological  Changes 

The  snake  pedal  model,  embedded  in  a  level-set  framework,  naturally  handles 
the  topological  changes  while  recovering  the  shapes  from  image  data.  We  introduce 
the  novel  concept  of  a  global  model  into  the  now  popular  PDE-based  curve/surface 
evolution  framework,  which  is  very  useful  in  the  shape  recognition  and  image  indexing 
tasks. 
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9.2    Future  Directions 

Our  experience  with  the  snake  pedal  model  has  revealed  several  issues  that 
are  relevant  to  the  continued  development  and  success  of  the  model.  This  section 
summarizes  these  issues  and  indicates  some  potential  further  research  directions. 

9.2.1  Invariant  Representation 

A  rigid  motion  invariant  representation  of  the  global  parameters  as  well  as  of  the 
coarse  level  coefficients  of  local  parameters  is  required  to  collect  the  statistics  using  a 
Bayesian  framework.  If  we  solve  the  global  parameters  in  accordance  with  a  certain 
predefined  order  while  fitting  the  model,  the  invariance  requirement  for  the  global 
parameters  can  be  satisfied.  Invariant  representation  of  the  local  representation  is 
currently  under  investigation  and  the  preliminary  results  are  encouraging. 

9.2.2  Introduce  Dvnamics  in  the  Representation 

Dynamics  is  necessary  in  the  motion  tracking  and  motion  synthesis  tasks.  Our 
physics-based  control  framework  can  be  extended  to  shape  and  non-rigid  motion 
estimators  by  incorporating  the  time  factor  into  the  governing  equation.  The  equation 
of  motion  will  become  a  second-order  PDE  and  efficient  numerical  techniques  need 
to  be  developed  to  solve  the  problem. 

9.2.3  Extend  the  hvbrid  geometric  active  model  to  3D 

Currently,  we  only  implemented  the  2D  hybrid  geometric  active  contours.  We 
need  to  extend  the  snake  pedal  model  to  cope  with  topological  changes  in  shape  in  3D 
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by  introducing  the  3D  hybrid  geometric  active  surfaces.  We  also  need  to  experiment 
with  more  image  data  in  3D. 


APPENDIX  A 
DERIVATION  OF  STIFFNESS  MATRIX  K 

In  Chapter  2,  we  obtained  the  internal  deformation  energy  in  an  active/deformable 
model.  To  obtain  the  expression  for  the  stiffness  matrix  K,  we  employ  the  finite  ele- 
ment method  (FEM)  to  discretize  the  model  in  material  coordinates  u  =  {u,v).  As 
shown  in  Fig.  5.2,  we  use  bilinear  quadrilateral  elements  in  conjunction  with  linear 
triangular  elements.  The  later  i  sused  for  representing  the  regions  around  the  north 
and  south  poles.  We  will  first  discuss  the  bilinear  quadrilateral  elements,  and  the 
linear  triangular  elements  can  be  handled  in  the  similar  way. 

A.l    Bilinear  Quadrilateral  Elements 

Bilinear  quadrilateral  elements  are  quite  common  in  FEM  literature  [87].  In 
this  section,  we  present  the  basic  transformation  equations  between  material  and 
element  coordinates.  Denote  the  finite  element  nodal  shape  functions  by  ,  where 
i  =  1, 2, no,  and  no  is  the  number  of  nodes  associated  with  each  element  Ej.  For 
the  bilinear  quadrilateral  elements,  we  have 

Ni{u,^)  =  ^{l+cocj,){l+4^iPi),  (A.l) 

where  {uji,ipi)  are  the  reference  coordinates  of  node  i  as  shown  in  Fig.  (A.l). 
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Figure  A.l.  Bilinear  quadrilateral  element,  where  {u,  v)  are  material  coordinates  and 
{lo,  ip)  are  the  reference  coordinates. 

The  reference  coordinates  (a;,  ^p)  are  related  to  the  material  coordinates  u  — 

{u,  v)  by 

2  2 

u  =  -{u-Uc)  iJ;  =  -{v-Vc),  (A.2) 

where  [uc,  Vc)  are  the  coordinates  of  the  element  center,  a  and  b  are  the  length  and 
width  of  the  element.  In  our  implementation,  we  use  a  =  b.  Differentiating  Eq.  (A.l) 
with  respect  to  u  and  v  gives 


(A.3) 


dNi     dNi  dto     dNi  dip      1  ^ 

-d^  =  a7    +  -di^d^  =  2^^^  + 


(A.4) 


.  ^  126 


When  integrating  function  f{u,  v)  over  element  Ej,  we  perform  the  integration  in  the 
reference  system: 

j      f{u,v)dudv  =  j  J  f{uj,ip)^dudip  (A.5) 


A.2    Derivation  of  the  Local  Stiffness  Matrix 

As  described  in  section  5.2.4,  we  use  the  membrane  plus  thin  plate  energy  ex- 
pressed in  Eq.(5.17)  as  the  internal  deformation  energy  imposed  on  the  snake, 

By  using  the  procedure  described  by  Terzopoulos  and  Metaxas  [78],  we  can  obtain 
the  local  stiffness  matrix.  For  the  membrane  energy  part,  if  we  denote  the  stiffness 
matrix  for  the  element  Ej  by       G  sjjnoxno^  ^-j^g  entries  of       can  be  represented  as 

{kDim  =     J  +  -^^)  du  dv,  (A.7) 

where  l,m  =  1,2,  ...,no.  Since  we  use  a  quadrilateral  element,  no  =  4  in  this  case. 
Using  the  formula  in  Eq.  (A.3),  (A. 4)  and  (A.5),  we  obtain  the  following  local  stiffness 
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matrix  for  the  membrane  energy 


D  a 


-2    -1  1 


-2 


1  -1 


-1  1 


1    -1  -2 


+  (--) 


1  -1 


-2  -1 


-1  -2 


(A.8) 


Adopting  the  similar  scheme,  we  obtain  the  following  stiffness  matrix  Kj  for  the  thin 
plate  energy  part  in  element  Ej 


*  ab 


1-1  1-1 


■1  1 


1-1  1-1 


-1      1-1  1 


(A.9) 


The  summation  of       and      gives  the  full  stiffness  matrix  for  element  Ej 


K^■  =  K4  +  K^ 


(A.IO) 


Note  that  both  and  K(  are  symmetric  positive  definite  (SPD)  matrixes,  therefore 
their  summation,  K-',  is  also  a  SPD  matrix. 

A. 3    Kronecker  Tensor  Product  Representation  of  the  Assembled  Stiffness  Matrix 

The  matrix  discussed  earlier  is  only  the  local  stiffness  matrix  for  element  Ej, 
which  describes  the  stress-strain  relationship  between  all  the  nodes  within  element 


128 


u 


1  2  3  4 


El 

5 

E2 

6 

E3 
7  1 

E4 

9 

Es 

10 

E6 
11  1 

V 

(a) 

Figure  A. 2.  A  small  3x4  mesh  in  {u,v)  space. 

Ej.  For  example,  ^3  =  4  means  applying  one  unit  displayment  to  node  3  results  in  4 
unit  force  on  node  2.  We  need  to  assemble  the  local  matrix  over  all  the  elements  to 
obtain  the  entire  stiffness  matrix  K.  To  achieve  this,  we  introduce  a  novel  method 
via  the  Kronecker  tensor  product  operation.  A  small  3x4  mesh  in  {u,  v)  space  as 
shown  in  Fig.  A. 2  is  used  to  illustrate  the  idea.  We  will  only  consider  the  first  part 
of  the  (the  first  term  in  the  summation  of  Eq.  (A. 8))  for  explanation,  the  second 
part  of  K4  and      can  be  handled  in  the  same  way. 

In  Fig.  (A. 2),  we  assign  a  number  to  each  node,  ranging  from  1  to  12,  and 
denote  each  element  by  Ei,  E2, E^.  We  choose  node  6  for  our  illustration,  and 
investigate  the  stress-strain  realtionship  between  node  6  and  every  other  node  in  the 
entire  mesh.  From  Eq.  (A. 8),  observe  that  applying  one  unit  displacement  to  node  6 
will  result  in  ^(^)  x  1  unit,  (^f )  J  x  (-1)  unit,  ^(  J)  x  (-2)  unit  force  on  node  2,  node 
1,  and  node  5,  respectively.  This  concept  can  be  visualized  by  the  nodal  molecule 
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Figure  A. 3.  Membrane  molecules 

representation,  as  shown  in  Fig.  (A. 3)  (a).  We  can  use  a  similar  representation  for 
elements  E2,  E4  and  £'5,  resulting  in  the  nodal  molecules  as  shown  in  Figs.  (A. 3)  (b), 
(c),  and  (d),  respectively.  Note  that  there  is  no  interaction  between  node  6  and  nodes 
in  elements  E3  and  Eq.  From  the  nodal  molecular  representation,  when  applying  one 
unit  displacement  to  node  6,  we  can  obtain  the  force  resulting  on  every  other  node  in 
the  entire  3x4  mesh  by  summing  up  all  the  molecular  values  at  the  same  molecular 
location.  The  result  is  shown  in  Fig.  (A. 4). 

The  assembled  molecular  representation  can  be  expressed  as  a  Kronecker  prod- 
uct fomula 
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Figure  A. 4.  Assembled  membrane  molecular  representation. 

The  assembled  stiffness  matrix  that  corresponds  to  the  first  part  of  the  mem- 
brane energy  (first  term  in  the  summation  of  Eq.  (A. 8))  is 


Kr"'  =  (y^)-A.®B 


where 


(A.12) 
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(A.13) 


We  now  prove  the  above  conclusion  by  expanding  the  Kronecker  tensor  product 
expression  in  (A.12).  According  to  the  Kronecker  product  definition  in  Eq.  (5.25), 
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In  this  matrix,  each  entry  represents  the  stress-strain  relation  between 

node  /  and  node  m,  where  l,m  =  1,2,  ...,12.  For  example,  the  entries  in  the  6th 
column  (A;^"'""^)(6  represent  the  force  caused  in  node  /  by  the  unit  displacement  in 
node  6.  We  can  see  from  column  6  that  one  unit  displacement  in  node  6  contributes 
(^|)  X  8  unit  force  to  itself,  and  (^J)  x  (-4)  unit  force  to  node  5  and  7,  (^^)  x  2 
unit  force  to  node  2  and  10,  (^^)  x  (-1)  unit  force  to  node  1,  3,  5  and  7,  and  zero 
force  to  node  4,  8  and  12.  This  is  exactly  what  Fig.  (A. 4)  conveys. 

For  the  second  part  of  the  membrane  energy,  performing  a  similar  procedure, 
we  have 

Kr"^==(^^).B.®A„  (A.15) 


where 
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Since  a  =  6  in  our  implementation,  the  entire  stiffness  matrix  corresponding  to  the 
membrane  energy  is  ,  • 


(A.17) 
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Similarly,  we  have  the  stiffness  matrix  corresponding  to  the  thin  plate  energy 


K,  =  (^)-A,®A,. 


(A.18) 


In  general,  for  an  m  x  n  mesh,  denote 


2  1 


1   4  1 


1   4  1 
1  2 


(A.19) 


Ay  — 


4  1 

1   4  1 


Bx  = 


1 

1  -1 
■1      2  -1 


1   4  1 
1  4 


e  3?" 


■1      2  -1 

-1  1 


(A.20) 


(A.21) 


134 


By  = 
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then  the  assembled  stiffness  matrix  corresponding  to  the  energy  expressed  in  Eq.(5.17) 


is: 


K  =       +  Kt  =  wi{A^  (g)  By  +      ®  Ay)  +  W2{A^  <8>  Ay), 


(A.23) 


where  u;i  and  W2  are  tension  and  rigidity  parameters.  Note  that  the  difference  be- 
tween the  first  column/row  and  last  column/row  in  and  A^  or  B^;  and  B^  is  due 
to  the  different  boundary  conditions,  namely,  the  natural  boundary  and  the  periodic 
boundary  conditions  respectively.  When  employing  the  same  boundary  conditions 
and  let  m  =  n,  A  =  Aj;  =  Aj,,  B  =  B^  =  B^,  Eq.  (A.23)  becomes 


K  =  w;i(A  (g)  B  +  B  O  A)  +  u;2(A  (g)  A). 


(A.24) 


This  is  the  same  equation  as  (5.23). 


APPENDIX  B 
RELATIONSHIP  OF     ~  </>i 


From  the  differential  geometry  theory,  we  have  the  relationship  between  the 
curvature  of  the  snake  pedal  kg  and  the  corresponding  higher  dimensional  function 


k.  =  V 


IIV02 


By  using  V02  =  JbV^i,  kg  can  be  expressed  as: 
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We  denote  each  entry  of  Jb  as: 


Jr  = 
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where  Js^^  =  J^^i  >  Js  is  the  symmetric  matrix.  Denote 
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Jb  can  be  written  as: 


JbV0i  = 
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The  magnitude  of  is: 


(B.6) 

Therefore,  we  can  represent  A;e  in  terms  of  by: 
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